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Abstract  of  Dissertation  Presented  to  the  Graduate  School  of 
the  University  of  Florida  in  Partial  Fulfillment  of  the  Requirements 
for  the  Degree  of  Doctor  of  Philosophy 

DYNAMICS  OF  INTERPLANETARY  DUST 
IN  THE  F CORONA 

by 

Edwin  T.  Rusk 
May  1986 

Chairman:  Nebil  Y.  Mi  scon i 

Major  Department:  Astronomy 

We  studied  the  dynamical  mechanisms  in  interplanetary  space  and 
in  the  F corona  by  numerical  simulations.  An  expression  for  the 
radiation  pressure  force  due  to  a rotating  spherical  source  of  radia- 
tion was  derived.  Also,  expressions  relating  the  variation  in 
inclination  and  the  longitude  of  the  ascending  node  to  the  solar 
magnetic  field  were  derived.  These  expressions  are  based  on  the 
spherical  source  surface  model  of  the  solar  magnetic  field. 

Simulations  of  particles  released  during  perihelion  passages  of 
comet  Encke  show  that  cometary  particles  have  lifetimes  SHORTER  than 
the  lifetime  calculated  by  Wyatt  and  Whipple  in  1950.  These  simula- 
tions also  resulted  in  higher  eccentricities  and  a definite  alignment 
of  the  particles'  aphelia  toward  a direction  20°  east  of  the  vernal 
equinox. 

An  expression  relating  the  size  of  a planet's  zone  of  influence 
to  perturbations  on  particles  in  solar  orbits  based  on  the  closest 
approach  between  the  planet  and  the  particle  show  that  the  expression 


v 


for  the  size  of  a planet's  zone  of  influence  is  not  singular,  but 
varies  with  the  particular  orbital  element  which  is  being  studied. 

Simulations  of  the  interaction  of  the  Lorentz  force  with  dust  in  the  F 
corona,  based  on  observed  solar  magnetic  field  values,  result  in  a 
SPREADING  OF  THE  INCLINATIONS  of  particles  in  circumsolar  orbits. 

This  result,  along  with  a reevaluation  of  recent  observations  of  the 
F corona,  leads  us  to  the  conclusion  that  the  shape  of  the  circumsolar 
dust  cloud  can  not  be  a ring,  but  must  be  a wide  band  or  a spherical 
shel 1 . 


CHAPTER  I 


INTRODUCTION 

On  a clear  dark  night  just  after  astronomical  twilight,  one  can 
observe  a section  of  the  sky  in  the  west  which  appears  brighter  than 
the  surrounding  regions.  This  light  comes  from  sunlight  reflected  by 
a myriad  of  small  particles,  less  than  .1  mm  across,  each  following 
its  own  orbital  path  about  the  sun.  Because  the  greatest  concentration 
of  these  particles  orbit  in  or  near  the  plane  of  the  earth's  orbit, 
and  thus  among  the  constellations  of  the  zodiac,  this  phenomenon  is 
called  the  zodiacal  light  and  the  interplanetary  particles  as  a group 
are  called  the  zodiacal  dust  cloud. 

The  F corona  is  the  innermost  region  of  the  zodiacal  cloud, 
culminating  in  a "dust  free"  zone  surrounding  the  sun.  The  light 
scattered  by  dust  particles  near  the  sun  forms  the  F corona  and  is 
part  of  the  zodiacal  light.  This  region  can  be  distinguished  from  the 
inner  K corona  by  the  spectra  of  the  light  emitted  by  these  two 
regions,  since  the  F corona  displays  the  Fraunhofer  spectrum  of  the 
sun  while  the  spectrum  of  the  K corona  is  characterized  by  the  emis- 
sion lines  typical  of  a very  high  temperature  plasma.  The  F corona 
is  the  final  resting  place  of  the  dust  grains  which  spiral  slowly 
inward  through  the  interplanetary  medium  due  to  the  radiation  induced 
Poynting-Robertson  drag  (henceforth  P-R  drag;  Poynting,  1903; 
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Robertson,  1937).  Upon  reaching  the  F corona,  the  particles  begin  to 
sublimate,  and  the  sublimation  rate  increases  until,  finally,  much 
smaller  remnants  of  the  original  dust  particles,  called  beta-meteoroids, 
are  blown  out  of  the  solar  system  by  radiation  pressure.  Because  the 
conditions  in  the  F corona  control  the  final  fate  of  interplanetary 
dust,  it  becomes  an  important  factor  in  determining  the  dust  budget, 
the  rate  at  which  zodiacal  dust  is  created  and  destroyed.  Unfortunately, 
at  the  present  time,  little  is  known  about  the  nature,  size,  composi- 
tion, or  dynamics  of  the  dust  particles  in  this  region. 

At  this  time,  the  source(s)  of  the  interplanetary  dust  are  not 
known.  Several  theories  have  been  put  forth;  the  most  popular  theories 
claim  that  the  dust  is  released  by  either  long  or  short  period  comets 
or  from  asteroid  and/or  meteoroid  collisions  (Gr’un  et  al.,  1985). 
Whichever  theory  is  true,  it  is  fairly  well  established  that  the  dust 
in  the  F corona  originated  farther  out  in  the  solar  system  (cf.  Mukai, 
1983)  and  reached  the  near  solar  regions  due  to  the  slow  decrease  in 
semimajor  axis  caused  by  the  P-R  drag. 

Scattering  measurements  (Greenberg  and  Gustafson,  1981)  indicate 
that  those  zodiacal  dust  particles  which  are  the  major  source  of  the 
zodiacal  light  are  "fluffy"  silicates,  irregular  in  shape  with  a density 
of  less  than  1 gm/cm  due  to  their  porous  nature.  Silicates,  and 
particularly  SiC^,  are  likely  candidates  for  F coronal  particles  since 
they  have  high  vaporization  temperatures  and  thus  can  survive  in  the 
coronal  regions.  Others  have  used  graphite  in  their  models  of  F 
coronal  dynamics  because  it  also  has  a high  vaporization  temperature, 
although  there  is  no  evidence  that  graphite  grains  comprise  a signifi- 
cant percentage  of  zodiacal  particles. 
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There  are  many  mechanisms  which  influence  the  motion  of  inter- 
planetary grains  as  they  orbit  the  sun.  Solar  radiation  acts  in  two 
ways  upon  orbiting  dust  particles.  First,  because  the  sun's  radiation 
field  is  radial  and  drops  off  as  the  inverse  square  of  the  distance 
from  the  sun,  the  force  generated  from  the  direct  pressure  of  solar 
radiation  forms  a constant  ratio  to  the  solar  gravitational  force,  and 
thus  the  combination  of  gravity  and  radiation  pressure  from  the  sun 
acts  upon  small  particles  exactly  as  if  the  sun  were  a non-radiating 
body  of  lesser  mass.  Secondly,  the  motion  of  each  particle  creates  a 
second  order  (in  v/c)  force  from  the  solar  radiation  due  to  both  the 
Doppler  shift  of  the  radiation  and  an  aberration  effect.  These  two 
second  order  effects  combine  to  form  the  P-R  drag.  This  drag  force 
opposes  the  motion  of  the  dust  particle,  reducing  the  particle's 
orbital  velocity  and  hence  the  kinetic  energy  of  its  orbit.  Since  the 
total  energy  of  an  orbiting  body  is  directly  related  to  its  semimajor 
axis,  one  effect  of  P-R  drag  is  to  slowly  decrease  the  semimajor  axis 
of  orbiting  dust  grains. 

Gravitational  effects  of  the  planets  cause  significant  perturba- 
tions upon  dust  particle  orbits.  Becuase  the  P-R  drag  continually 
reduces  the  semimajor  axis,  an  orbiting  particle  will  have  its  orbital 
path  crossing  those  of  each  of  the  inner  planets  at  some  point  during 
its  lifetime.  This  increases  the  likelihood  of  one  or  more  close 
planetary  encounters.  The  nature  of  the  effects  of  gravitational 
interactions  with  the  planets  is  therefore  different  for  dust  par- 
ticles as  opposed  to  large  bodies  such  as  meteoroids  or  asteroids, 
which,  if  they  are  in  orbits  which  to  not  intersect  those  of  the 
planets,  will  not  have  these  "catastrophic"  planetary  encounters. 
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The  interplanetary  magnetic  field  is  an  extension  of  the  solar 
magnetic  field  into  the  planetary  environment  carried  by  the  solar  wind. 
It  will  alter  the  orbital  path  of  small  charged  dust  particles  by  the 
Lorentz  force.  The  interplanetary  magnetic  field  can  be  considered 
to  be  a combination  of  two  components,  a radial  component  directed 
inward  to  or  outward  from  the  sun  and  a tangential  component  perpen- 
dicular to  both  the  radial  direction  and  the  direction  of  the  solar 
rotation  axis. 

In  the  region  of  the  F corona,  the  radial  component  of  the  inter- 
planetary magnetic  field  is  much  stronger  than  the  tangential  component. 
Here,  the  magnetic  field  becomes  an  important  factor  in  perturbing  the 
particle's  orbit  because  the  usually  stronger  P-R  drag  is  balanced  by 
a loss  of  particle  size  due  to  sublimation  which  results  from  the 
high  temperatures  of  this  region.  This  balance  allows  the  weaker 
Lorentz  force  to  play  a more  significant  role.  Outside  the  F corona, 
the  Lorentz  force  becomes  negligible  compared  to  P-R  drag.  Because 
the  radial  component  of  the  magnetic  field  drops  off  as  the  square  of 
the  distance  from  the  sun,  as  does  P-R  drag,  the  Lorentz  force  due 
to  the  radial  component  of  the  interplanetary  field  can  be  neglected 
everywhere  except  in  the  F corona.  The  tangential  component  of  the 
interplanetary  magnetic  field  only  decreases  as  the  heliocentric 
distance  to  the  1st  power,  so  at  large  distances  from  the  sun,  the 
Lorentz  force  due  to  this  tangential  component  is  of  the  same  order 
as  the  P-R  drag  for  micron-sized  particles.  This  is  of  significant 
importance  in  determining  the  influx  of  interstellar  dust  into  the 
solar  system  (Gustafson  and  Misconi,  1979). 
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In  the  F corona,  the  radiation  acceleration  to  gravitational 
acceleration  ratio  will  no  longer  remain  constant  when  the  particles 
reach  the  region  where  they  begin  sublimation.  The  magnitude  of  the 
acceleration  due  to  radiation  pressure  increases  as  the  particle  size 
decreases  and  results  in  the  particle  experiencing  a continual  decrease 
in  the  effective  mass  of  the  sun.  This  equates  to  an  increase  in  the 
potential  energy  of  the  grain,  which  results  in  an  increase  in  the 
semimajor  axis  of  the  particle's  orbit  (Lamy,  1974;  Mukai , 1983).  We 
will  call  this  the  sublimation  effect.  At  some  specific  distance 
from  the  sun,  the  sublimation  rate  is  such  that  it  will  balance  the 
P-R  drag  and  the  semimajor  axis  of  a particle's  orbit  will  temporarily 
remain  stable.  This  region  should  be  observable  as  an  enhancement  in 
the  brightness  of  zodiacal  light  due  to  the  enhancement  of  particle 
number  density. 

Recent  observations  have  detected  such  an  enhancement  in  the 
density  of  orbiting  dust  particles  at  about  4 RQ(sol ar  radii).  Maiharaetal. 
(1985)  used  photoelectric  scans  to  detect  thermal  emission  from  dust 
in  the  F corona  during  the  June  1983  total  solar  eclipse  in  Indonesia. 

They  described  their  findings  as  dust  rings  around  the  sun  at  3.5  RQ 
and  4 RQ  near  the  solar  equatorial  plane.  During  the  eclipse  of 
November  1966  Peterson  (1967,  1969)  observed  an  emission  peak  in  the 
infrared  [X  = 2.2  ym)  surrounding  the  sun  at  about  4 RQ.  During  the 
same  eclipse,  MacQueen  (1968)  observed  brightness  enhancements  at  3.4 
Rq  and  4.1  RQ  in  a radial  scan,  also  at  2.2  pm,  and  at  4 Rq  during  a 
balloon  flight  on  January  9,  1967.  Finally,  during  the  February  1979 
eclipse  Beavers  et  al . (1980)  detected  a radial  velocity  component  of 


-6- 


the  dust  in  the  F corona  which  matches  that  expected  from  Keplerian 
orbital  motion. 

Because  the  sun  is  not  a point  source,  the  radiation  from  different 
points  of  origin  on  the  solar  surface  will  appear  at  different  fre- 
quencies due  to  the  Doppler  effect.  The  P-R  drag  and  radiation  pressure 
comprise  the  non-Doppler  average  of  these  source  points,  and  when  they 
are  subtracted  from  the  total  radiation  distribution,  the  residual 
radiant  energy  will  add  an  additional  drag  term  to  the  P-R  drag  (Guess, 
1962).  This  is  because  the  leading  hemisphere  of  the  solar  disk  appears 
bluer,  and  hence  the  particle  intercepts  photons  of  higher  energy  from 
the  leading  hemisphere  as  opposed  to  the  trailing  hemisphere. 

Another  effect  of  the  sun's  finite  diameter  is  due  to  the  Doppler 
shift  of  solar  radiation  caused  by  the  sun's  rotation.  This  effect, 
similar  to  the  effect  of  a finite  solar  diameter,  will  generate  a 
force  in  the  solar  equatorial  plane  normal  to  the  line  between  the 
particle  and  the  sun's  center.  Unless  the  dust  particle  is  orbiting 
in  the  solar  equatorial  plane,  this  force  will  have  a component  per- 
pendicular to  the  particle's  orbital  plane  and  can  thus  affect  the 
inclination  of  the  particle's  orbit. 


CHAPTER  II 


PLANETARY  PERTURBATIONS  AND  THE  STRUCTURE  OF 
THE  ZODIACAL  CLOUD 

A.  Introduction 

Cometary  debris  and  particles  released  in  asteroidal  collisions 
are  believed  to  be  the  two  main  sources  of  the  zodiacal  dust.  Because 
the  asteroids  are,  for  the  most  part,  in  relatively  circular  orbits 
at  low  inclinations  to  the  invariable  plane  of  the  solar  system,  it 
is  expected  that  these  particles  would  form  a homogeneous,  time  in- 
variant component  of  the  zodiacal  cloud.  Much  work  has  been  done 
(cf.  Gustafson,  1985)  on  the  dynamics  of  dust  grains  in  circular 
orbits  under  the  influence  of  planetary  perturbations  combined  with 
P-R  drag.  Although  certain  features,  such  as  the  Schuerman  arcs 
(Schuerman,  1980)  which  result  from  orbital  resonance  with  the 
planets,  are  expected  to  form  regions  of  enhanced  dust  density,  in 
general,  asteroidal  dust  is  expected  to  manifest  itself  as  a fairly 
homogeneous  dust  cloud. 

Comets,  on  the  other  hand,  have  highly  eccentric  orbits  which 
frequently  undergo  drastic  changes  due  to  close  interactions  with  the 
planets,  especially  Jupiter.  Most  of  the  dust  from  comets  is  released 
while  the  comet  is  near  perihelion,  where  the  higher  temperature 
allows  the  ices  which  bind  the  refractory  material  to  the  comet  to 
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subl i mate , thus  releasing  this  material.  A change  in  the  comet's 
perihelion  distance  will  therefore  strongly  affect  the  amount  of  dust 
released.  If  the  zodiacal  dust  is  mainly  due  to  cometary  debris,  it 
may  be  possible  to  determine  which  comet  released  the  dust  and  when 
it  was  released,  if  structures  relating  to  specific  perihelion  pas- 
sages of  likely  cometary  suppliers  of  the  interplanetary  dust  can  be 
observed. 


B.  The  Method 

To  determine  the  effect  of  gravitational  perturbations  on  the 
cometary  dust,  196  particles  were  released  at  differing  perihelion 
passages  of  comet  Encke.  Because  the  orbital  period  of  this  comet 
has  a seven  to  two  commensurabi 1 ity  with  the  orbital  period  of  the 
planet  Jupiter,  particles  were  chosen  in  groups  of  seven  consecutive 
perihelion  passages  in  order  to  determine  if  the  evolution  of  the 
cometary  particle  orbits  depend  on  the  Jupiter-sun-particle  angle 
at  release.  Each  of  these  seven  release  "phases"  could  then  be 
studied  separately. 

Due  to  the  commensurabi 1 i ty  in  orbital  period,  comet  Encke' s 
orbital  parameters  are  strongly  influenced  by  the  planet  Jupiter. 

The  result  of  this  influence  is  a periodic  variation  in  the  comet's 
orbital  elements  with  respect  to  Jupiter's  orbital  plane  (Whipple, 
1940).  Figure  2.1  (Whipple,  1940,  p.  730)  demonstrates  how  these 
parameters  vary  with  time.  In  our  simulations,  we  used  the  periodic 
variation  in  inclination  to  determine  the  times  at  which  our  particles 
would  be  released.  By  studying  particles  released  when  the  comet  was 
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Variation  of  the  orbital  parameters  of  comet  Encke.  Here,  the  orbital 
parameters  are  measured  with  respect  to  the  orbital  plane  of  Jupiter. 
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at  a minimum  possible  inclination,  at  a midpoint  (10°)  of  inclination, 
and  at  a maximum  inclination,  we  could  see  if  there  were  any  effects  of 
initial  inclination  on  the  final  fate  of  each  cometary  particle.  Also, 
it  was  important  to  include  a range  of  input  values  to  get  an 
accurate  determination  of  the  distribution  of  cometary  dust  at  any 
given  time. 

At  release,  we  assumed  that  the  relative  velocity  between  the 
comet  and  the  particle  was  zero.  Any  micron-sized  particle  leaving 
a comet  is  suddenly  subjected  to  a radiation  pressure  induced  accelera- 
tion much  greater  than  that  on  the  comet,  and  its  velocity  (the  comet's 
velocity)  is  too  high  to  remain  on  the  comet's  orbital  path.  Because 
radiation  pressure  effectively  decreases  the  central  force,  the  dust 
particle's  orbit  increases  in  both  semimajor  axis  and  eccentricity. 

In  fact,  if  the  radiation  acceleration  to  gravitational  acceleration 
ratio,  8,  is  such  that 


where  ec  is  the  eccentricity  of  the  comet,  then  the  dust  particle  will 
be  released  in  a hyperbolic  orbit.  This  sets  a lower  limit  to  the  size 
of  cometary  dust.  In  our  simulations,  we  chose  to  use  dust  particles 
30  microns  in  diameter  (8  = 0.02)  to  represent  an  intermediate  value 
for  zodiacal  dust,  which  is  believed  to  be  between  1 and  100  microns 
(Giese  and  Grun,  1976)  in  size.  The  resulting  particles'  orbital  paths 
had  semimajor  axes  of  2.89  A.U.  and  eccentricities  of  .884,  compared 
to  comet  Encke's  values  of  2.21  A.U.  and  .848,  respectively.  The 
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other  orbital  parameters  remain  as  those  of  the  comet,  except  for  the 
orbital  period  which  increased  by  nearly  50%.  The  significance  of  the 
increase  in  semimajor  axis  is  that  the  particle's  orbit  crosses  that 
of  Jupiter,  so  that  close  encounters  with  this  massive  planet  become 
likely.  This  can  greatly  affect  the  subsequent  evolution  of  the 
particle's  orbit. 

From  the  equations  of  Wyatt  and  Whipple  (1950),  the  time  it  would 
take  a particle  to  reach  the  inner  solar  system  from  the  initial  input 
semimajor  axis  and  eccentricity  can  be  determined.  This  has  been 
tabulated  for  several  values  of  the  particle's  semimajor  axis  (table 
2.1).  It  can  be  seen  that  a particle  requires  between  15,000  and 

20.000  years  to  reach  the  inner  solar  system,  and  so  the  times  of 
particle  release  in  our  simulations  were  set  between  26,000  BC  and 

6.000  BC  (figure  2.2).  The  longer  timespan  of  20,000  years,  as  com- 
pared to  the  5,000  years  required  to  cover  all  possible  release  times 
for  particles  solely  under  the  influence  of  P-R  drag,  was  needed  to 
compensate  for  the  effect  of  planetary  perturbations  on  the  spiral  ling- 
in  time  of  the  cometary  dust. 


Table  2.1.  The  time  required  for  a particle  released  from  comet  Encke 
to  reach  a particular  semimajor  axis  under  the  influence 
of  the  Poynting-Robertson  effect,  and  the  resulting 
eccentricity. 


Semi major 
Axis  (AU) 

Eccentrici ty 

Time  from 
release  (yr) 

1.050000 

0.705010 

13082.38 

0.950000 

0.678761 

14008.29 

0.850000 

0.647662 

14968.27 

0.750000 

0.610368 

15964.00 

0.650000 

0.565083 

16994.10 

0.550000 

0.509462 

18052.78 

0.450000 

0.440659 

19124.28 

0.350000 

0.355914 

20172.94 
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comet.  Points  not  on  the  curve,  resulting  from  an  error  in  the  input  computa- 
tions, will  not  apply  to  P.  Encke,  but  their  inclusion  does  not  affect  the 
conclusions  drawn  in  this  chapter.  The  curve  represents  the  inclination  of 
comet  Encke  during  this  period. 
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The  evolution  of  the  particle's  path  through  the  solar  system  was 
solved  by  using  Encke's  method.  Here  the  osculating  orbit  of  the 
particle  was  divided  into  two  hundred  equal  time  increments,  and,  at 
each  of  these  times,  the  acceleration  of  the  dust  particle  from  P-R 
drag  as  well  as  that  due  to  the  planets  from  Venus  through  Jupiter  was 
determined.  The  perturbations  on  the  osculating  orbit  were  computed 
using  the  Gaussian  perturbation  equations  (Moulton,  1914).  The 
changes  in  the  orbital  elements  were  then  integrated  over  one  orbit 
using  Simpson's  rule,  and  the  osculating  orbit  was  corrected  at  the 
end  of  each  orbit.  When  tests  were  performed  on  the  method  using 
particles  run  without  planetary  encounters,  the  constant  resulting 
from  the  effect  of  radiation  pressure  as  determined  by  Wyatt  and 
Whipple  (1950), 


C = 


ae  -V5  (1  . 2, 

P P P 


5 


remained  constant  to  at  least  four  decimal  places  (the  extent  of  our 
printout)  thoughout  the  20,000  orbits  of  the  test  particle. 

When  the  particle  approached  within  the  zone  of  influence  of  a 
planet  (Chapter  IV  will  contain  more  about  the  range  of  planetary 
influence)  the  integration  was  completed  and  the  particle's  path  was 
computed  by  determining  the  change  in  position  and  velocity  vectors 
in  Cartesian  coordinates  (Cowell's  method).  The  size  of  the  time 
increment  was  decreased  to  1/40, 000th  of  an  orbital  period  throughout 
the  Cowell's  method  procedure.  This  method  was  continued  until  the 
particle  was  sufficiently  far  from  the  planet  that  Encke's  method 
could  once  again  be  used. 
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Each  particle  was  then  run  until  either  the  semimajor  axis  fell 
below  .3  astronomical  units,  the  perihelion  distance  fell  below  two 
solar  radii,  or  the  particle's  semimajor  axis  failed  to  fall  below  2.5 
A.U.  after  10,000  orbits,  which  indicated  that  the  particle's  orbit 
became  dominated  by  Jupiter.  Most  of  the  particles  required  between 
10,000  and  20,000  orbits  to  reach  one  of  these  criteria  and  took 
approximately  one  hour  of  cpu  time  per  1000  orbits;  therefore  each 
particle  run  required  about  one  day  to  complete  on  a VAX/750. 

C.  Results 

The  distribution  of  the  particles  was  analyzed  by  generating 
histograms  of  several  orbital  elements  as  well  as  the  time  since  re- 
lease for  values  of  the  semimajor  axis  from  2.75  to  .35  A.U.  or  for 
20  different  dates  separated  by  1500  year  intervals.  Also,  the 
inclination  of  each  particle  was  plotted  as  a function  of  its  longi- 
tude of  the  ascending  node.  These  figures  can  be  found  in  Appendix  B. 

In  the  plots  of  inclination  vs.  longitude  of  the  ascending  node 
there  appears  to  be  no  correlation  between  the  two  orbital  elements. 

In  fact,  the  orbital  planes  of  the  dust  particles  appear  to  be  dis- 
tributed randomly  about  the  ecliptic  (i  = 0°)  plane,  with  the  number 
density  of  particles  falling  off  at  inclinations  greater  than  about 
20°  (figure  2.3).  However,  a glance  at  the  histograms  of  inclination 
show  a marked  decrease  in  the  particle  density  for  inclinations 
between  0°  to  3°.  This  is  due  to  the  nature  of  measuring  an  angle 
such  as  inclination,  which  is  restricted  to  the  range  from  0°  to  180°. 
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particle  is  represented  by  up  to  25  points  corresponding  to  the  25 
values  of  the  semimajor  axis  studied. 
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A histogram  across  all  semimajor  axes  (figure  2.4)  strongly  exhibits 
this  effect. 

Statistically,  if  a particle  is  placed  in  a random  orbit,  the 
angular  momentum  vector  of  a particle's  orbit  should  be  equally  likely 
to  point  in  any  direction.  Specifically,  this  means  that  the  proba- 
bility of  the  angular  momentum  unit  vector  ending  in  any  surface 
element  of  a unit  sphere  is  the  area  of  that  surface  element  divided 
by  the  area  of  the  sphere  (4n).  If  a sphere  is  sliced  parallel  to  its 
equatorial  plane  such  that  the  areas  of  the  resulting  surfaces  are 
equal,  the  size  of  the  angle  subtended  by  the  areas  is  not  constant. 

In  order  to  have  a constant  area  on  the  sphere  the  slices  must  be  made 
at  equal  increments  of  the  cosine  of  the  inclination.  Histograms  of 
the  number  density  of  particles  at  equal  intervals  of  the  cosine  of 
the  inclination  for  various  values  of  the  semimajor  axis  (Appendix  B, 
figures  B.26  to  B.50)  show  that  the  particles  are  indeed  driven 
toward  lower  inclinations. 

The  number  densities  of  the  particles  found  at  the  lowest  inclina- 
tions (figure  2.5)  imply  that  the  particles  are  selectively  removed 
from  those  ranges  of  inclination  inhabited  by  the  planets.  This  is 
to  be  expected  since  close  encounters  between  the  dust  particles  and 
the  planets  are  much  more  likely  to  occur  when  the  orbital  planes  of 
the  two  bodies  nearly  coincide.  In  figure  2.5,  the  three  outermost 
planets  in  our  study,  Jupiter,  Mars,  and  Earth,  are  found  in  the 
zones  of  the  last  two  bars,  and  the  cosine  of  the  inclination  of  the 
planet  Venus  is  .99824,  which  puts  that  planet  squarely  in  the  fourth 
bar  from  the  right,  a bar  that  is  clearly  low  in  particle  number 
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figures  B.l  through  B.25. 
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inclinations  near  the  ecliptic. 
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density.  This  close-up  study  of  the  inclination  does  not  detract  from 
the  fact  that,  overall,  particle  density  increases  toward  lower  values 
of  the  inclination.  Unfortunately,  there  is  not  enough  structure  in 
the  available  data  to  determine  whether  the  symmetry  plane  of  the 
particles  in  the  simulations  is  related  to  the  orbital  plane  of  any 
particular  planet. 

It  has  been  frequently  stated  that  particles  in  the  inner  solar 
system  will  be  in  nearly  circular  orbits  due  to  the  decrease  in 
eccentricity  resulting  from  the  P-R  drag,  and  Wyatt  and  Whipple  (1950) 
have  been  referenced  as  the  source  of  this  belief.  However,  table  2.1 
demonstrates  that  cometary  particles  will  be  found  in  eccentric  orbits 
even  when  their  semimajor  axes  are  as  small  as  0.35  A.U.  Figure  2.6 
shows  how  P-R  drag  alone  would  affect  the  evolution  of  both  the  peri- 
helion and  aphelion  distances  of  a particle  released  from  comet  Encke. 
Observe  that  the  distance  between  the  two  values  stays  large  for  most 
of  the  particle's  lifetime. 

The  eccentricity  of  the  cometary  simulation  particles  has  been 
catalogued  in  Appendix  B,  for  values  of  the  semimajor  axis  from 
2.75  A.U.  to  0.35  A.U.  in  steps  of  .1  A.U.  Figure  2.7  shows  the  median 
eccentricity  as  a function  of  semimajor  axis  for  our  196  particles. 

The  lower  curve  in  this  figure  represents  the  single  value  which  would 
result  if  the  effect  of  the  planets  were  excluded.  BY  CONTRAST y IT 
CAN  THUS  BE  SEEN  THAT  THE  PLANETS  INCREASE  THE  ECCENTRICITY  OF  THE 
MAJORITY  OF  THE  DUST  PARTICLES  RELEASED  FROM  COMET  ENCKE. 

The  effect  of  the  planets  on  the  time  the  particles  required  to 
fall  toward  the  sun  can  be  seen  in  figure  2.8.  The  upper  curve  in 
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Figure  2.7.  The  median  eccentricity  of  the  simulation  particles  as  a function  of 

semimajor  axis.  The  lower  curve  represents  the  effect  of  P-R  drag  alone. 
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this  figure  represents  the  infall  time  derived  from  P-R  drag.  IT  IS 
CLEAR  THEN  THAT  THE  MAJORITY  OF  THE  PARTICLES  HAVE  THEIR  LIFETIMES 
SIGNIFICANTLY  SHORTENED  BY  INTERACTIONS  WITH  THE  PLANETS.  This  result 
does  not  show  whether  planetary  encounters,  as  well  as  slow  secular 
changes  due  to  more  distant  planets,  affect  the  particle's  semimajor 
axes  directly  to  cause  this  reduction  in  lifetime  or  indirectly 
through  a strengthening  of  the  P-R  drag  due  to  the  increase  in  eccen- 
tricity discussed  above. 

The  P-R  lifetime  of  a particle  in  an  eccentric  orbit  is  smaller 
than  the  lifetime  of  a particle  with  the  same  semimajor  axis  in  a 
circular  orbit  (Wyatt  and  Whipple,  1950),  because  the  increase  in  the 
P-R  drag  near  perihelion  more  than  compensates  for  the  decrease  when 
the  particle  is  at  distances  greater  than  its  semimajor  axis.  A par- 
ticle that  has  had  its  eccentricity  increased  by  planetary  encounters 
will  have  its  P-R  lifetime  correspondingly  decreased;  thus  the  effect 
of  the  planets  on  the  eccentricities  of  the  particles  can,  at  least 
partly,  explain  the  reduction  in  P-R  lifetime. 

To  determine  whether  this  effect  is  the  sole  cause  of  the  de- 
crease in  the  lifetime  of  the  simulation  particles,  the  eccentricity 
and  "age"  of  the  particles  were  compared  with  the  P-R  values  based  on 
the  median  eccentricity  of  a particle  at  2.45  A.U.  (figures  2.9  and 
2.10)  and  at  1.55  A.U.  (figures  2.11  and  2.12).  The  semimajor  axis 
2.45  A.U.  was  chosen  as  the  largest  value  of  semimajor  axis  that 
guaranteed  that  the  particles  would  not  cross  Jupiter's  orbital  path. 
The  other  value,  1.55  A.U.,  was  chosen  since  a small  deviation  in 
the  smoothness  of  the  upper  curve  in  figures  2.7  and  2.9  occurred 
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of  the  upper  curve  of  2.45  A.U. 
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Figure  2.10.  Same  as  figure  2.8  with  the  exception  noted  in  figure  2.9. 
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Figure  2.11.  Same  as  figure  2.9  except  the  lower  curve  starts  at  1.55  A.U. 
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Figure  2.12.  Same  as  figure  2.10  with  the  lower  curve  starting  at  1.55  A.U.,  as  in 
figure  2.11. 
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between  1.55  and  1.65  A.U.  The  resulting  figures  (2.10,  2.12)  show 
that  even  after  the  higher  eccentricity  is  compensated  for,  the 
particles  still  show  a tendency  to  reduce  their  semimajor  axes  more 
rapidly  than  P-R  drag  would  predict.  THIS  IMPLIES  THAT  THE  effect  of 
GRAVITATIONAL  INTERACTIONS  WITH  THE  PLANETS  MOST  OFTEN  TEND  TO  REDUCE 
A PARTICLE'S  SEMIMAJOR  AXIS. 

The  longitudes  of  the  perihelia  of  the  particles  in  our  simula- 
tions are  preferentially  aligned  away  from  the  direction  <L  = 240°, 
measured  in  a prograde  sense  from  the  direction  of  the  vernal  equinox. 
This  can  be  seen  in  the  histograms  of  Appendix  B (figures  B.101  through 
B.125)  and  has  been  summed  over  all  semimajor  axes  in  figure  2.13. 

THIS  ALIGNMENT  IS  NOT  RELATED  TO  THAT  OF  THE  COMET , WHICH  HAS  ITS  VALUE 
OF  THE  LONGITUDE  OF  PERIHELION  VARYING  THROUGH  AN  ENTIRE  CYCLE  OF  360° 
DURING  THE  20,000  YEAR  RELEASE  PERIOD.  Also,  the  perihelion  directions 
appear  to  be  more  strongly  ordered  for  a constant  value  of  the  semi- 
major axis  rather  than  across  all  semimajor  axes  at  any  given  time, 
implying  that  such  an  alignment  would  be  a stable  feature  in  three 
dimensions. 


D.  Conclusion 

Any  results  based  on  the  dust  released  from  comet  Encke  should 
apply  to  any  prograde  comet  of  low  inclination,  since  it  is  apparent 
from  the  simulations  that  the  planet  Jupiter  causes  a catastrophic 
realignment  of  the  dust  when  the  dust  is  in  Jupiter  crossing  orbits. 
Even  those  dust  particles  with  high  inclinations  have  a good  chance 
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of  a close  encounter  with  Jupiter,  since  P-R  drag  will  eventually 
bring  the  particles'  orbital  paths  to  intersection  with  Jupiter's 
orbit. 

One  effect  of  planetary  perturbations  is  to  slowly  draw  the  dust 
particles  toward  the  planets'  orbital  planes.  This  is  in  line  with 
observations  (Leinert  et  al . , 1981)  which  show  the  very  well  known  fact 
that  zodiacal  light  brightness  increases  toward  the  ecliptic. 

A surprising  result  is  the  determination  that  the  eccentricities 
of  the  particle  orbits  remain  high  throughout  much  the  particles' 
lifetimes.  Observations  have  not  yet  been  made  which  can  distinguish 
between  particles  in  circular  orbits  and  those  in  eccentric  orbits. 

The  most  recent  results  (East  and  Reay,  1984)  confirm  that  at  least 
50X  of  the  dust  particles  are  in  closed  prograde  orbits  about  the  sun. 
Theoretical  studies  of  the  dynamics  of  particles  near  the  sun  (Lamy, 
1974;  Mukai  et  al . , 1974)  imply  that  most  interplanetary  grains 
eventually  lose  enough  mass  through  sublimation  to  evolve  into  par- 
ticles with  hyperbolic  orbits  (3-meteoroids),  thus  explaining  the 
percentage  which  are  not  in  closed  orbits. 

The  asymmetry  in  the  longitude  of  perihelion  of  the  cometary 
dust  is  apparently  a feature  of  the  solar  system  in  general  and  not 
specific  to  comet  Encke,  since  the  orientation  of  the  comet  appears 
to  have  no  effect  on  the  direction  of  the  feature.  That  there  has 
been  no  observation  of  such  a feature  would  imply  that  the  majority 
of  interplanetary  dust  particles  are  in  orbits  of  low  eccentricity. 

It  appears  that  dust  from  comet  Encke  does  not  meet  all  the 
criteria  required  of  a source  of  the  zodiacal  cloud.  The  median 
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lifetime  of  particles  released  from  the  comet  is  shortened,  so  that 
the  total  particle  mass  required  to  supply  the  zodiacal  dust  would  be 
even  more  than  previous  estimates  based  solely  on  the  P-R  lifetime 
(Mukai  et  al . , 1983),  which  conclude  that  cometary  dust  is  an  insuf- 
ficient source.  Radiation  pressure  eliminates  cometary  dust  below  the 
lower  limit  set  by  equation  2.1,  thus  much  of  the  mass  of  a periodic 
comet  will  be  ejected  from  the  solar  system,  further  reducing  the 
initial  supply  of  interplanetary  dust.  Finally,  the  lack  of  asymmetry 
in  the  distribution  of  interplanetary  dust  implies  that  the  dust  orbits 
are  predominantly  of  low  eccentricity,  a result  we  have  shown  to  be 
inconsistent  with  a cometary  origin. 


CHAPTER  III 


INITIAL  PARTICLE  DISTRIBUTION  OF  THE  F CORONA 

In  order  to  study  particle  dynamics  in  the  F corona,  the  first 
order  of  business  is  to  determine  the  size,  charge,  chemical  composi- 
tion, and  orbital  elements  of  particles  as  they  enter  the  region.  We 
can  consider  the  F corona  and  the  rest  of  the  zodiacal  cloud  as  two 
separate  regions  if  our  criterion  for  dividing  the  two  is  to  consider 
the  rest  of  the  zodiacal  cloud  as  that  region  where  radiation  pressure 
(including  P-R  drag)  and  gravitational  interactions  are  the  only 
significant  forces  affecting  the  dynamics  of  interplanetary  dust. 

The  Lorentz  force  from  the  interplanetary  magnetic  field  cannot 
be  entirely  neglected  anywhere  in  the  solar  system,  but  for  10-100  pm 
particles  it  is  small  compared  to  the  perturbations  from  the  planets. 

In  particular,  in  the  outer  regions  of  the  solar  system  where  the 
tangential  component  of  the  interplanetary  field  generates  forces 
comparable  to  the  P-R  drag,  the  effect  of  the  planets  is  signifi- 
cantly enhanced  due  to  their  large  masses  and  the  weakening  of  the 
sun's  gravitational  field. 

The  size  of  dust  particles  in  the  F corona  depends  on  several 
factors  which  can  reduce  the  size  of  particles  as  they  approach  the 
sun,  as  well  as,  of  course,  the  initial  sizes  of  the  interplanetary 
particles.  Particle  collection  studies  in  the  Earth's  upper  atmosphere 
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(Brownlee  et  al.,  1980)  as  well  as  interplanetary  space  probes  (Pailer  and 
Griin,  1980)  and  the  study  of  lunar  microcraters  (LeSargent  and  Lamy, 

1978)  point  to  there  being  two  populations  of  particles  responsible 
for  the  zodiacal  light.  Population  I particles  are  dark  fluffy  con- 
glomerates of  silica  and  silicates  with  most  of  the  particles  being 
between  10  and  100  microns  in  diameter.  These  particles  are  in  pro- 
grade, low  eccentricity  orbits  about  the  sun.  The  particles  which 
comprise  population  II  are  less  than  1 micron  in  diameter  and  appear 
to  be  streaming  outward  from  the  sun  in  hyperbolic  orbits.  It  is  be- 
lieved that  the  number  density  of  particles  in  the  two  populations  are 
approximately  equal  (Giese  and  Gr'un,  1976). 

Apparently,  the  population  II  particles  are  the  remnants  of  the 
larger  particles  after  the  dust  has  undergone  sublimation  in  the  F 
corona.  This  would  seem  to  be  the  most  reasonable  explanation, 
although  this  has  not  been  discussed  in  the  literature.  It  is  based 
on  this  conjecture  that  we  intend  to  study  only  those  particles  larger 
than  1 micron.  In  any  case,  it  has  been  shown  that  population  I par- 
ticles are  responsible  for  90%  of  the  zodiacal  light  (East  and  Reay, 
1984). 

The  reduction  in  size  of  the  interplanetary  particles  can  result 
from  four  effects: 

1.  Sublimation  of  the  particle, 

2.  Sputtering,  which  is  an  eroding  of  the  particle  caused  by 
bombardment  of  solar  wind  ions, 

3.  Rotational  bursting,  and 

4.  Collisions. 
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Of  these  four  mechanisms,  sublimation  can  be  ignored  since  this 
mechanism  is  confined  (by  definition)  to  the  F corona.  In  the  case  of 
rotational  bursting,  there  are  many  mechanisms  known  which  result  in 
the  rotation  of  dust  particles  (Radzievski i , 1954;  Misconi,  1975;  Paddack 
and  Rhee,  1975) . Because  the  particles  appear  to  be  fluffy  conglomerates 
and  not  solid  crystals,  their  tensile  strength  is  relatively  low  and 
thus  the  particles  are  apparently  easy  to  burst.  Misconi  (1975)  has 
shown  that  random  changes  in  rotation  due  to  solar  wind  ion  impacts 
may  cause  rotational  bursting  of  particles  during  solar  flares, 
although  mechanisms  which  slow  the  particle  spin  rate  could  not  be 
included  in  the  analysis.  Other  rotation  mechanisms  are  capable  of 
spinning  the  particles  up  to  bursting  speed  much  more  rapidly  than 
random  ion  impacts,  but  these  depend  on  the  existence  of  specific  types 
of  structure  on  the  surfaces  of  the  dust  particles,  as  well  as  a stable 
alignment  of  the  particle  spin  axis.  Such  mechanisms  are  not  well 
understood  at  the  present  time. 

In  our  present  work,  we  will  not  include  rotational  bursting  as 
a significant  factor  in  determining  the  size  of  the  dust  grains  as 
they  reach  the  F corona.  There  are  three  main  reasons  for  excluding 
rotational  bursting: 

1.  A search  for  enhancement  of  the  zodiacal  light  from  solar 
flares  failed  to  discover  any  evidence  of  such  an  occurrence. 

2.  A particle  breaking  up  will  experience  a sudden  change  in 
the  central  force  field  due  to  an  enhancement  in  radiation 
pressure,  thus  many  such  particles  will  be  driven  from  the 
solar  system  (see  equation  2.1). 

3.  Misconi's  (1975)  analysis  implies  that  bursting  will  occur 
only  in  particles  less  than  1 micron  in  diameter  based  on 
a tensile  strength  for  solid  materials  (for  interplanetary 
grains  the  value  may  be  much  lower).  Applicability  of  this 
mechanism  to  the  larger  grains  in  our  study  is  unclear. 
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Sputtering,  like  sublimation,  causes  a slow  degradation  in  the 
sizes  of  interplanetary  dust  particles.  Unlike  sublimation,  sputter- 
ing can  affect  particles  outside  the  F corona.  The  sputtering  rate  of 
interplanetary  dust  particles  measured  in  A/year  (Bandermann,  1968)  is 

ds  = 3.2  2A 

pa 

3 

where  dt  is  measured  in  years,  p is  measured  in  grams/cm  , and  a is 
measured  in  A.U.  Here  we  have  taken  a value  intermediate  to  Bandermann' s 
maximum  and  minimum  values  for  the  sputtering  rate  of  iron  particles  and 
adjusted  for  the  lower  density  of  fluffy  particles.  An  adjustment  was 
also  made  in  the  sputtering  yield  per  ion  as  being  a factor  of  the 
density  of  iron.  We  can  then  calculate  the  size  of  a particle  in  a 
CIRCULAR  orbit  by  substituting  for  dt  from  the  equation 

da  = (t)  dt  3.2 

a 

o c 

(Wyatt  and  Whipple,  1950),  where  a = -4-.  Here  p is  measured  in  grams 

bp 

per  cm3,  and  s is  in  Angstroms,  as  was  ds.  Combining  equations  3.1 
and  3.2  we  get 

ds  _ 8 da 
s’"?  a ’ 

which,  when  integrated  gives 
a 8/9 

s = S0  (f-)  . 3.3 

0 

If,  for  example,  sQ  = 30  y at  1 A.U.,  then  it  will  have  been  reduced 
to  2.1  y at  .05  A.U.  (approximately  10  RQ). 
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Collisions  will  certainly  reduce  the  size  of  interplanetary 
particles.  In  our  study  of  the  F corona  we  will  ignore  the  effect 
of  collisional  reduction  since  particles  broken  up  by  collision  will 
most  likely  become  3-meteoroids  and  never  reach  the  F corona. 

The  end  result  of  the  question  of  the  size  of  dust  particles 
entering  the  F corona  is  that  sputtering  is  the  main  reducing  factor 
of  particle  size  in  the  inner  solar  system.  Based  on  equation  3.3, 
we  expect  the  interplanetary  dust,  which  is  mainly  between  10  and  100 
microns  in  size  at  1 A.U.,  to  be  between  .7  and  7 microns  in  size 
upon  reaching  the  F corona. 

There  is  much  disagreement  as  to  the  charge  of  particles  in  the 
F corona.  Burke  and  Silk  (1974)  suggest  a negative  charge  for  graphite 

particles  at  temperatures  less  than  8000  K.  The  applicability  of  this 

result  must  be  questioned  because  of  the  low  photoelectric  yield  in 
the  environment  under  study  (HI  regions)  as  compared  to  the  near  solar 
environment.  Millet  et  al . (1980)  derive  a value  of  - +3  volts  for  the 
potential  of  spherical  graphite  particles  near  the  sun,  after  reaching 
a negative  value  of  -40  volts  at  about  10  RQ.  Mukai  (1981)  states 
that  the  particle  will  remain  positively  charged  throughout  the  inner 
solar  system,  with  the  potential  gradually  increasing  above  +10  volts 
in  the  minimum  solar  wind  while  reaching  only  +4.8  volts  if  the  solar 
wind  is  at  maximum.  All  of  these  results  are  for  graphite,  although 
Burke  and  Silk  state  that  the  behavior  of  silicates  is  closely  related. 

Lacking  a definitive  value  of  the  charge  of  dust  grains  in  the  F 

corona,  we  will  assign  an  arbitrary  value  of  +10  volts.  The  implica- 
tions of  an  error  in  our  chosen  value  of  the  charge  will  be  discussed 
in  the  next  chapter. 
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Si nee  the  chemical  composition  of  the  dust  particles  has  been  deter- 
mined (Brownlee  et  al. ,1980),  at  least  for  those  particles  in  the 
neighborhood  of  the  Earth,  as  being  silica  and  silicates,  there  remains 
only  to  determine  the  orbital  elements  of  the  particles  as  they  enter 
the  F corona. 

In  the  previous  chapter,  the  orientation  of  the  particles'  orbital 
planes  with  respect  to  the  ecliptic  plane  was  determined  to  vary,  with 
somewhat  more  particles  having  low  inclination  orbits,  rather  than 
high  inclination  orbits.  There  appears  to  be  no  preferential  value 
for  the  longitude  of  the  ascending  node  of  the  particles  studied.  We  will 
adopt  Lei nert  et  al. 1 s (1976)  fan  model  for  the  particle  density  since 
it  fits  our  results  rather  nicely. 

The  semimajor  axis  of  F coronal  dust  will  be  adopted  as  being 
between  3.5  RQ  and  4.5  RQ,  from  the  observations  discussed  in  Chapter 
I and  from  theoretical  considerations  (Lamy,  1974). 

The  eccentricities  of  dust  particle  orbits  will  be  reduced  by  P-R 
drag  as  given  by  Wyatt  and  Whipple  (1950).  Planetary  effects  on  the 
eccentricity  of  the  orbits  are  smal 1 compared  to  that  of  P-R  drag  in- 
side the  orbit  of  Mercury.  Table  3.1  is  an  extension  of  table  2.1  and 
shows  the  eccentricity  of  an  average  cometary  particle  as  it  approaches 
the  F corona.  It  can  be  seen  from  table  3.1  that  the  eccentricity  of 
a particle  in  the  F corona  is  less  than  0.04,  even  for  our  cometary 
particle  with  a high  eccentricity  at  release.  In  our  study,  we  will 
assume  that  our  particles  are  initially  in  circular  orbits,  and  so 
there  is  no  need  to  consider  the  asymmetry  in  the  direction  of 
perihel ion. 
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Table  3.1.  The  effect  of  P-R  drag  on  the  eccentricity  and  the  P-R 
lifetime  of  a cometary  particle  near  the  sun. 


Semimajor 
Axis  (A.U.) 

Eccentricity 

Time  from 
Release  (yrs.) 

0.350000 

0.355914 

20172.94 

0.350000 

0.307176 

20667.81 

0.250000 

0.254588 

21127.42 

0.200000 

0.199120 

21535.42 

0.150000 

0.142485 

21874.59 

0.100000 

0.087226 

22128.76 

0.080000 

0.066262 

22203.70 

0.060000 

0.046378 

22262.47 

0.040000 

0.027986 

22304.68 

0.030000 

0.019543 

22319.50 

0.025000 

0.015563 

22325.32 

0.020000 

0.011776 

22330.09 

CHAPTER  IV 


PLANETARY  INFLUENCE  IN  THE  F CORONA 

Besides  determining  the  distribution  of  dust  particles  at  the 
beginning  of  their  "lifetime"  in  the  F corona,  we  need  to  know  just 
what  effect  planetary  perturbations  have  on  particles  oribiting  in 
the  F corona.  Since  planets  are  in  nearly  circular  orbits  and  the 
region  of  the  F corona  is  within  10  solar  radii  (.05  A.U.)  of  the 
center  of  the  sun,  a planet's  semimajor  axis  is  the  approximate  dis- 
tance between  the  planet  and  the  particle  at  all  times.  Thus  by 
determining  the  size  of  the  region  in  which  a planet  can  have  a sig- 
nificant influence  on  an  orbiting  particle,  we  will  determine  if  any 
such  planet  has  a zone  of  influence  as  large  as  its  semimajor  axis. 

Tisserand's  (1889)  definition  of  the  "sphere  of  action"  of  a 
planet  is  based  on  the  equality  of  tidal  vs.  gravitational  accelera- 
tion ratios  of  the  sun  and  planet.  Op i k (1951)  and  others  based  their 
relation  on  equating  the  differential  solar  and  planetary  forces  on 
the  particle.  Neither  expression  was  formulated  to  describe  the  zone 
of  influence  surrounding  a planet  when  considering  the  small,  but 
significant,  perturbative  effects  of  the  planets  on  a particle's 
orbital  elements.  For  the  purpose  of  determining  these  effects  on 
interplanetary  dust  we  derive  a zone  of  influence  based  on  equating 
the  gravitational  forces  of  the  sun  and  planet  and  demonstrate  its 
applicability. 
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In  the  process  of  computing  the  orbital  evolution  of  dust  par- 
ticles in  interplanetary  space,  it  was  necessary  to  determine  whether 
or  not  a particle  was  close  enough  to  a planet  that  integrating  the 
perturbing  forces  using  Encke's  method  was  no  longer  valid.  Both 
expressions  for  the  "sphere"  of  influence  of  a planet  proved  to  be 
inadequate  when  the  geometries  of  subsequent  interactions  were  con- 
sidered. In  particular,  when  each  relation  was  scaled  to  one  planet, 
it  could  not  be  generalized  to  the  other  planets.  Consequently,  the 
formula  for  the  zone  of  influence  must  depend  strongly  on  the  specific 
subject  to  be  studied.  We  have  pursued  this  problem  from  several  dif- 
ferent viewpoints  and  expect  that,  besides  being  useful  in  studies  of 
interplanetary  dust,  our  findings  may  apply  to  studies  of  asteroids, 
comets,  and  the  stability  of  asteroidal  satellites. 

The  general  expression  for  the  radial  distance  to  the  boundary  of 
the  planetary  zone  of  influence  on  small  particles  can  be  written  as 

s = k-f(a  )*g(m  , other  planetary  parameters)  4.1 

r r 

where  k is  a constant  which  depends  on  the  nature  of  the  problem  under 
study,  f and  g are  functions  which  may  also  have  the  same  dependence 
as  k,  ap  is  the  planet's  semimajor  axis,  and  nip  is  the  planet's  mass. 
The  other  parameters  may  include  the  eccentricity  of  the  planet,  as 
well  as  the  other  orbital  parameters  and  such  factors  as  oblateness 
and  axial  tilt.  Functions  f and  g are  generally  not  separable  but  we 
have  found  that  the  effect  of  semimajor  axis  can  be  separated  from  the 
other  parameters.  The  equations  of  Ti sserand  (1889)  and  Opik  (1951) 


can  be  written  as 
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1 2/5 

S = (M1)  apmp  7 5 and 

0 K v 

1 1 1/3  1/3 

s = 2 ^2M~^  apmp  » respectively. 


4.2 

4.3 


If  we  equate  the  gravitational  and  radiation  forces  of  the  sun  and 
planet,  a different  relationship  is  obtained.  Consider  a particle  of 
mass  m at  a distance  s from  a planet  of  mass  nip  and  semimajor  axis  ap, 
and  at  a distance  r from  the  sun.  The  magnitude  of  the  force  on  the 
particle  due  to  the  planet  is 


Gm  m 

F =-E 


2 * 


4.4 


while  the  magnitude  of  the  force  on  the  particle  due  to  the  sun  is 


F = GM  m ^ , 

0 r2 


4.5 


where  3 is  the  ratio  of  radiation  pressure  to  the  sun's  gravitational 

force.  Since  ap  - s < r < ap  + s,  we  can  write  r as  r = ap  + 6s  where 

6 is  a number  between  -1  and  +1  which  depends  on  the  planet-particle-sun 

angle.  Substituting  = Mq(1  - 3)  we  can  set  the  ratio  of  the  solar 

2 

force  to  the  planetary  gravitational  force  equal  to  some  constant  k , 
where 


k2Gm„m 


GM'm 

0 


(ap  + 6s ) 


2 , which  reduces  to 


m_  1/2  m 1/2  , 

- kan(-E)  [1  - (-£)  r1  . 

P m'  M‘ 


4.6 


s 


4.7 
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Equation  4.7  can  be  approximated  by  a simple  power  law  when  the 
mass  of  the  planet  and  k are  small  (compared  to  the  sun),  giving 


M 


■Wp1'2 


4.8 


Compare  this  with  the  results  of  Tisserand  and  Opik,  equations  4.2 
and  4.3. 

A computer  simulation  technique  was  used  to  find  the  most  applicable 

expression  for  our  study.  Each  encounter  was  run  using  a 10,000  step 

Cowell's  method  computation  over  one  synodic  period  of  the  orbit. 

For  a particle  in  circular  orbit  at  10  A.U.,  the  errors  in  eccentricity 

-5 

and  semimajor  axis  incurred  over  one  orbit  were  less  than  1 x 10  , and 

-9 

5 x 10  A.U.,  respectively,  and  there  was  no  observable  error  in 

inclination.  The  particle  was  initially  in  a circular  orbit  and  the 

planet  was  given  a fixed  orbital  eccentricity.  The  encounter  occurred 

at  perihelion  or  aphelion,  depending  on  whether  the  particle's  orbit 

was  inside  or  outside  the  planet's  orbit.  The  particle's  orbital  plane 

coincided  with  that  of  the  planet  except  when  the  change  in  inclination 

of  the  particle  was  under  study.  In  this  case,  the  particle's  initial 

inclination  was  determined  by  the  distance  of  closest  approach,  so 

that  encounters  above  and  below  the  planet's  orbital  plane  could  be 

studied.  A comparison  was  made  between  the  minimum  encounter  distance 

and  a or  m under  the  condition  of  constant  Ae,  Ai , or 
P P a 

Figure  4.1a  shows  the  effect  of  a„  on  the  zone  size  for  a con- 
3 P 

stant  Ae  and  verifies  the  linear  relation,  common  to  all  derivations, 
between  the  zone  boundary  and  the  semimajor  axis.  This  result  also 
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shows  that  the  uncoupling  of  the  functional  forms  of  ap  and  mp  in 
equation  4.1  is  valid.  The  constant  change  Ae  = 0.0014  chosen  in 
deriving  figures  4.1a  and  4.1b  (and  later  4.3)  is  the  minimum  value 
that  guarantees  that  the  zones  of  influence  of  the  inner  five  planets 
do  not  overlap.  Similar  linear  relationships  (not  shown)  were  found 


when  Ai  and  were  held  constant. 

a 

The  encounters  were  computed  with  the  gravitational  interactions 
of  the  sun  and  planet  affecting  the  particle  continuously  throughout 
one  synchronous  orbit  and  not  just  at  the  point  of  closest  approach  to 
the  planet.  Equations  4.2,  4.3,  4.7,  and  4.8  were  derived  by  consider- 
ing only  one  point  in  the  particle-planet-sun  system.  So  far,  we  have 
assumed  a relation  between  the  distance  of  closest  approach  in  our 
simulation  results  and  the  distance  s in  these  equations.  The  linear 
relationship  between  ap  and  the  distance  of  closest  approach  in  figure 
4.1a  demonstrates  that  our  assumption  is  useful.  Some  effects  of  the 
difference  between  the  derived  equations  and  the  numerical  simulations 
of  the  complete  encounters  will  be  discussed  later. 

Figure  4.1b  shows  the  effect  of  the  planet's  mass  on  the  zone 
boundary  when  Ae  is  held  constant.  The  general  trend  follows  the  power 
law  of  equation  4.8  and  diverges  from  the  power  law  in  the  same  manner 
as  equation  4.7.  Figure  4.2  is  a graphic  representation  of  equation 
4.7.  While  this  and  the  previous  figure  show  a striking  similarity, 
the  spread  due  to 
mD  1/2 

[1  - k6(-*)  ] 1 

m; 


is  greater  than  the  spread  observed  in  the  simulations.  This  is  a 


CLOSEST  APPROACH  DISTANCE  (IN  AU]  CLOSEST  APPROACH  DISTANCE  {IN  AU| 
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Figure  4.1.  Relation  between  closest  approach  distance  and 

planetary  characteristics  for  a constant  change  in 
a particle's  orbital  eccentricity,  a)  Relation 
between  closest  approach  distance  and  the  semimajor 
axis  of  a planet  of  1 earth  mass,  b)  Relation  be- 
tween closest  approach  distance  and  the  mass  of  a 
planet  orbiting  at  10  A.U. 
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Figure  4.2.  A graphic  representation  of  equation  4.7  for  a planet  orbiting  at  10  A.U. 
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consequence  of  the  difference  mentioned  above,  due  to  the  fact  that, 
during  an  encounter,  6 assumes  values  other  than  +1  or  -1.  However, 
the  power  law  relation,  found  by  averaging  the  two  curves,  is  not 
affected  by  this  consideration  and  demonstrates  the  applicability  of 
equations  4.7  and  4.8. 

Figure  4.3  shows  the  effect  of  planetary  eccentricity  (ep)  on  the 
zone  boundary.  As  implied  in  equation  4.1,  this  effect  cannot  be 
separated  from  the  functional  dependence  on  planetary  mass.  This  is 
evidenced  by  the  fact  that  these  results  match  those  for  high  mass 
planets  but  diverge  from  those  of  low  mass.  This  is  evidently  caused 
by  a difference  in  the  velocity  relationship  between  the  particle  and 
the  planet.  When  the  particle  encounters  an  exterior  planet  at  peri- 
helion, the  particle  will  usually  have  the  greater  velocity.  However, 
if  the  encounter  distance 


s < 


2ap6P 
Tl  + en) 


the  planet  will  have  the  greater  velocity  and  will  overtake  the  par- 
ticle. In  figure  4.3,  the  inequality  is  true  when  s is  less  than 
.95  A. U.  A similar  effect  applies  when  the  particle's  position  is  out- 
side that  of  a planet  at  aphelion.  The  resulting  increase  in  the 
encounter  time  thus  causes  the  downward  curvature  at  the  low  mass  end 
of  the  curves  in  figure  4.3. 

Figures  4.4a  and  4.4b  show  the  relationship  between  planetary  mass 
and  the  distance  of  closest  approach  under  the  condition  that  — = 

d 

5 x 10"5  (figure  4.4a)  or  Ai  = 0.0011  (figure  4.4b).  Unlike  the  result 
for  a constant  Ae  in  figure  4.2,  the  boundary  follows  the  form 
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Figure  4.3.  This  is  the  same  as  figure  4.1b  except  the  planet's  orbit  has  an 
eccentricity  of  0.05. 
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Figure  4.4.  Relation  between  the  mass  of  a planet  and  the 

closest  approach  distance  for  a constant  change 
in  a particle's  a)  semimajor  axis,  b)  inclination. 
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m 3/8 


4.9 


a result  which  matches  no  theoretical  formula  yet  derived.  Note,  also, 
the  divergence  from  the  power  law  for  planets  of  high  mass,  similar  to 
the  form  in  equation  4.7. 

THE  FORM  OF  THE  RELATION  DETERMINING  THE  ZONE  OF  INFLUENCE  TO 
HIGHLY  DEPENDENT  ON  THE  NATURE  OF  THE  PERTURBING  EFFECT  TO  BE  OB- 
SERVED. Note  that  the  expressions  most  applicable  to  interactions  with 
orbiting  particles  involve  powers  of  the  planetary  mass  not  formerly 
considered,  and  the  size  of  the  zone  boundary  depends  upon  the  planet's 
eccentricity.  The  importance  of  the  3/8  exponent  in  the  power  law 
relation,  discovered  for  constant  Ai  and  suggests  a need  for  further 

a 

analytical  study. 

We  can  gain  insight  into  the  significance  of  planetary  perturba- 
tion on  the  orbital  elements  of  the  particles  by  considering  the  "k" 
factor  in  the  expression  for  the  size  of  the  zone  boundary  as  a 
measure  of  the  inverse  magnitude  of  the  perturbations.  Substituting 
the  semimajor  axis  in  place  of  the  size  of  the  zone  boundary  in  equa- 
tion 4.7,  and  replacing  the  exponent  (1/2)  by  the  more  generalized  £, 


which  implies  that,  in  the  F corona,  the  relative  influence  of  the 
planets  depends  solely  on  the  mass  of  each  planet.  There  is  the  factor 
of  1 + 6 to  be  considered,  since  6 is  always  very  close  to  -1.  In  the 
F corona,  1 + 6 is  equal  to  the  angular  distance  between  the  sun  and 


1 mn  E 

Perturbations  = v = (-^)  (1  + 6)  ; 


4.10 
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the  particle,  as  seen  from  the  planet.  If  we  consider  the  maximum 

angular  distance  in  radians  between  the  sun  and  the  particle,  then  that 

angle  is  the  radius  of  the  particle's  orbit  divided  by  the  semimajor 

4R 

axis  of  the  planet's  orbit.  Entering  the  value  — - for  1 + 6,  we  arrive 

ap 

at 


m e 

4Rn 


4.11 


The  significance  of  this  result  is  that  perturbations  due  to  the 
nearest  planet,  Mercury,  are  less  than  those  of  Venus,  Earth,  and 
Jupiter,  in  spite  of  the  fact  that  these  planets  are  much  farther  from 
the  F corona.  In  fact,  JUPITER  IS  THE  STRONGEST  PERTURBER  OF  ORBITING 
DUST  IN  THE  F CORONA , REGARDLESS  OF  THE  EXPONENT  OF  THE  MSS.  From  the 
simulations  above,  we  know  that,  for  a constant  change  in  inclination, 
the  exponent  E in  equation  4.11  is  3/8.  In  this  case,  Mercury,  Venus, 
and  Earth  all  contribute  at  least  half  of  the  inclination  change  caused 
by  Jupiter.  Table  4.1  displays  the  relative  effects  of  the  planets 
(Jupiter  = 1)  at  the  maximum  angular  distance  between  the  sun  and  a 
particle  orbiting  at  4 solar  radii,  for  several  values  of  the  exponent 
of  the  mass,  E , in  equation  4.11. 


Table  4.1.  The  relative  magnitudes  (Jupiter  = 1)  of  planetary  perturbations  in  the  F corona,  for 
various  values  of  the  dependence  on  planetary  mass. 
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CHAPTER  V 


THE  EFFECT  OF  THE  SOLAR  MAGNETIC  FIELD 
A.  An  Overview 

The  solar  magnetic  field,  carried  by  the  solar  wind,  is  expected 
to  interact  with  orbiting  particles  through  the  Lorentz  force.  The 
most  recent  models  of  the  solar  magnetic  field  (Schulz  et  al . , 1978; 
Hoeksema,  1984)  assume  that  there  exists  a surface  above  the  solar 
photosphere  where  all  the  magnetic  field  lines  are  radial,  and  that 
beyond  this  surface  the  magnetic  field  streams  outward  at  the  velocity 
of  the  solar  wind,  describing  an  Archemedian  spiral  (Parker,  1963). 
Inside  this  "source  surface"  the  field  increases  in  complexity  until, 
at  the  photosphere,  it  can  be  described  by  a series  of  multipole 
components.  The  source  surface  gives  accurate  fits  when  it  is 
modelled  between  2 and  2.85  solar  radii.  Although  ellipsoidal  shapes 
have  been  used  for  the  source  surface  (Schulz  et  al.,  1978),  the  spherical 
model  of  Hoeksema  (1984)  used  in  this  analysis  is  equally  accurate  in 
reproducing  the  observations  and  can  be  used  to  derive  analytical 
expressions. 

Since  we  will  be  dealing  with  the  magnetic  field  outside  the 
source  surface,  we  must  apply  the  equations  for  the  Parker  spiral 
pattern.  The  radial  component  of  the  magnetic  field  at  a distance  r 
from  the  center  of  the  sun  is 
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5.1 


where  R$  is  the  radius  of  the  source  surface  and  Bq  is  the  magnetic 
field  strength  at  the  source  surface,  taking  the  finite  velocity  of 
the  field  into  account  when  determining  the  position  on  the  source 
surface  where  the  field  originated.  The  tangential  component  of  the 
field  at  the  heliocentric  distance  r is 


where  e is  the  colatitude  of  the  position  of  field  measurement.  The 

«0RS 

factor  — - can  be  calculated  as  follows.  Approximately,  the  sun 


5.2 


or 


5.3 


rotates  once  each  month  and  there  are  about  n x 10'  seconds  in  a 


year,  so  that  the  angular  rotation  rate  of  the  sun 


271  x -12  - = 2.4  x 10“6  rad/sec 
TT  X 10 


5 

and  Rs  = 2.5  RQ  where  RQ  is  about  7.5  x 10  km.  Thequiet  solar  wind 
velocity  is  400  km/sec  (Gustafson  and  Misconi,  1979)  so 


2.4  x 10 


4 x 10 


,-6 


x 7.5  x 10 


2 


5 


= 4.5  x 10“  3 


The  maximum  value  of  the  tangential  component  for  any  given  radial 
distance  to  the  sun  occurs  in  the  solar  equatorial  plane  where  sin  e = 1. 
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In  the  F corona,  the  distance  to  the  sun  is  less  than  10  solar  radii, 
so  the  maximum  value  of  the  tangential  component  of  the  interplanetary 
magnetic  field  is 

B (r)  < B Jr)  x 4.5  x 10'3  x 9 x 1 
$ r 

-2 

which  is  <5x10  times  the  field  strength  of  the  radial  component  and 
therefore  will  be  ignored  in  the  analytical  analyses  that  follow  in 
this  chapter. 

All  of  the  previous  expressions  require  knowing  the  field  strength 
at  the  source  surface.  Although  I have  ascribed  no  functional 
dependence  to  Bq,  it  is  in  fact  a function  of  time  and  the  position 
angles  e and  <j>.  The  magnetic  field  at  the  source  surface,  can  be 
found  from  the  multi  pole  components  by  substituting  r = Rs  into  the 
expression  for  the  radial  component  of  the  solar  magnetic  field  used 
by  Hoeksema  (1984,  p.  44,  equation  3.1)  in  determining  these  com- 
ponents. The  resulting  expression 
R 1+2 

Bq  = (2Z+1)(-^)  P^  (cos  e)  (gZmcos  rrn j>  + h^sin  m$)  5.4 

is  time  dependent  because  the  multi  pole  components,  g^  and  h^  are  them- 
selves time  dependent.  The  functions  are  the  second  order  Legendre 
functions,  the  first  three  terms  of  which  can  be  found  in  Appendix  A. 

The  Lorentz  force  acting  on  the  orbiting  dust  grains  due  to  the 
solar  magnetic  field  will  perturb  the  orbits  of  the  coronal  dust.  In 
general,  this  force  will  act  on  any  orbiting  body  and  the  force  acting 
on  that  body  for  a given  position,  velocity,  and  time  is 
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- V 


SW 


) X B 


5.5 


where  vg  is  the  velocity  of  the  body  (grain)  with  respect  to  the 
inertial  coordinate  system,  v$w  is  the  solar  wind  velocity  with  respect 
to  that  system,  q is  the  net  charge  of  the  orbiting  body,  $ is  the 
magnetic  field  vector  for  the  given  position  and  time,  and  c is  the 
velocity  of  light.  Note  that  vsw  also  represents  the  velocity  of  the 

magnetic  field  since  the  field  is  "frozen"  in  the  plasma. 

The  dust  particles  in  this  region  are  expected  to  be  in  nearly 
circular  orbits  due  to  the  decrease  in  eccentricity  from  the  P-R  drag 
(discussed  in  Chapter  III),  and  therefore  their  velocity  vectors  are 
nearly  perpendicular  to  the  magnetic  field  lines.  The  Lorentz  force 
is  then  mainly  due  to  the  velocity  of  the  particle,  since  the  solar 
wind  velocity  is  nearly  parallel  to  the  magnetic  field  in  this  region. 
The  skewing  of  the  field  lines  (i.e.,  the  Archemedian  spiral;  Parker, 
1963),  as  shown  from  the  above  analysis,  takes  place  at  further 
heliocentric  distances  than  those  under  consideration.  The  resulting 
force  is  normal  to  the  particle's  orbital  plane,  so  only  those  orbital 

elements  which  can  be  affected  by  a force  perpendicular  to  the 

particle's  orbit  will  display  significant  changes. 


B.  The  Dipole  Component 

During  sunspot  minimum,  the  solar  magnetic  field  can  be  charac- 
terized by  a polar  dipole  field.  The  higher  order  multipole  components 
are  very  weak  compared  to  the  dipole  terms,  and  the  magnitude  of  the 
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equatorial  dipole  component  of  the  field  is  low  enough  so  that  the 
total  dipole  field,  which  is  a sum  of  the  polar  and  equatorial  dipole 
components,  is  inclined  less  than  5°  to  the  solar  rotation  axis. 

For  approximately  half  of  a solar  cycle,  there  is  a strong  equa- 
torial dipole  term,  especially  during  those  times  around  sunspot 
maximum.  The  equatorial  dipole  is  usually  accompanied  by  relatively 
strong  components  in  the  higher  order  multipoles. 

The  radial  field  at  the  source  surface  due  only  to  these  dipole 
terms  can  be  written  (from  equation  5.4)  as 

R 3 1 

B0  = 3 (r^)  l P^(cose) (g1  cosmij)  + hlmsinm<t>)  , 5.6 

s ffpO 

where  e is  the  polar  angle  and  <j>  is  the  azimuthal  angle  as  measured 
from  the  direction  of  zero  solar  longitude.  Expanding  this  gives  only 
three  terms,  since  sin(0<j>)  = sin(O)  = 0.  Substituting  P^(cose  ) from 
Appendix  A,  equations  A. 9,  we  then  get 

R 3 

BQ  = 3 (^)  [g1Qcose  + gnsin  e cos  <}> 

+ h^sine  sin^]  . 5.7 

Since  the  tangential  component  of  the  interplanetary  magnetic 
field  can  be  neglected,  we  can  write  the  expression  for  the  magnetic 
field  at  the  particle's  position  by  substituting  equation  5.7  into 
equation  5.1 

Bf(r)  = k(r)  [g1Qcose  + gnsine  cos<}> 


+ h^sine  si n^]  , 


5.8 
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R R 2 

where  k(r)  = 3(-^-)(— ) . Using  the  transformations  of  Appendix  A 
s 

(equations  A. 13),  we  can  write  the  expression  for  the  magnetic  field 
in  terms  of  the  particle's  orbital  inclination  i,  the  longitude  of  the 
ascending  node  ft,  and  u,  the  position  angle  of  the  particle  from  its 
ascending  node.  Note  that  the  orbital  elements  of  the  particle  used 
here  and  in  future  expressions  are  based  on  the  solar  equatorial 
coordinate  system.  The  expression  for  the  magnetic  field  at  the 
particle  becomes 

Br(r)  = kMfg^gSinftsini  - g^cosftcosu  - sinftsinucosi) 

- h^(sinftcosu  + cosftsi nucosi )]  . 5.9 

The  force  on  the  particle  due  to  the  magnetic  field  can  be  found 
by  substituting  equation  5.9  into  the  Lorentz  force  equation  (5.5). 
Because  the  magnetic  field  is  assumed  to  be  radial,  expression  5.5  can 
be  rewritten  as 

qvQ 

W = Br(r)  5.10 

where  W is  the  component  of  the  force  perpendicular  to  the  particle's 
orbital  plane.  The  other  components  of  the  force  (S  and  T)  are,  in 
this  case,  0. 

To  determine  the  effect  of  this  force  on  the  particle's  orbital 
elements,  we  examine  the  Gaussian  perturbation  equations  (Moulton, 
1914;  see  Appendix  A,  equations  A.14-A.19).  It  is  apparent  that  only 
the  inclination,  ascending  node,  and  the  longitude  of  perihelion  can 
be  affected  by  the  Lorentz  force.  Since  the  eccentricity  is  assumed 


-58- 


to  be  zero,  there  is  no  perihelion  direction  and  we  need  only  concern 
ourselves  with  changes  in  i and  ft. 

In  order  to  better  understand  the  role  each  component  of  the 
magnetic  field  plays  in  altering  the  orbital  elements  of  a dust  particle 
in  the  F corona,  the  expressions  for  the  change  in  the  orbital  elements 
due  to  components  g^g,  g^,  and  h^  will  be  written  separately.  The 
effect  of  g^  on  the  inclination  is 


If  we  substitute  for  dt  in  terms  of  u and  du  (Appendix  A,  equation 
A. 13)  and  integrate  both  sides  over  one  orbit,  with  the  assumption 
that  the  variations  in  i and  ft  are  small  over  the  course  of  one  orbit, 
we  get 

| d i = -C(r)  g1Q  ^sinucosw  d u , 
and 

2 

$ dft  = -C(r)  g^g  $ sin  u d u , 


5.11a 


and  the  effect  of  g^g  on  the  ascending  node  is 


5.11b 


where  C(r)  = (-7 )k(r),  which  can  be  integrated  to  give 
n^acm 


Ai(gig)  0 


5.12a 


AQ(g10)  - -irC(r)  g1Q  . 


5.12b 
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The  effect  of  the  dipole  component  of  the  solar  magnetic 
field  on  a particle's  orbital  alignment 

. qv 

dl  = ^naan^1^  (cosficos^  ~ sinfisinucosucosf)g11  5.11c 

and 

dfi  qv 

dl  = ^nacm^^r^  (cos^sinwcosucscf  - sinasin  ucotf) 5. lid 

results  in  the  integral  equations 

2 

j di  = C(r)  g^  [cosft  ^ cos  u d u - sinftcosf  ^ sinucosu  du] 
and 

Q 

j dQ  = C(r)  g.j  [cosficscf  f sinucosu  du  - sinftcoti  f sin  u du] 
which,  after  performing  the  required  integrations,  give 

Ai(gn)  = ttC ( r ) gn  cos n 5.12c 

and 

Afi(g1]L)  = -TiC(r)  gn  sin^cotf  . 5.12d 

Finally,  we  take  the  expressions  for  the  change  in  the  orbital 
elements  due  to  the  h^  component, 

d . qv 

df  = (nacm^(r)  hn  (sinacos  m + cosGsinucosucosf ) 5. lie 

and 

, qv 

df  = (nacm^r)  hll  (sinfisinwc0s^csci  + cosGsin  ucoti)  5.  Ilf 


which  we  can  integrate  as  before  to  get 
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Ai(hn)  =irC(r)  hn  sinfi  5.12e 

and 

Aft(h^)  =irC(r)  h^  cosftcoti  . 5.12f 

From  equation  5.12a,  it  is  evident  that  the  polar  dipole  com- 
ponent (g^g)  has  no  effect  on  the  inclination  of  a particle's  orbital 
plane  and  from  equation  5.12b,  it  is  seen  that  there  is  a constant 
monotonic  change  in  ft  due  to  the  polar  dipole  component.  Because  the 
elements  i and  fi  are  based  on  a chosen  coordinate  system,  the  results 
for  the  effect  of  the  polar  dipole  component  can  be  extrapolated  to 
show  that  the  effect  of  any  solar  dipole  field  is  to  precess  the 
particle's  orbital  plane  about  the  dipole  axis.  When  one  or  both  of 
the  equatorial  dipole  terms  are  strong,  the  precession  results  in  a 
significant  change  in  the  inclination.  Keep  in  mind  that  this  is  the 
effect  of  a radial  field  which  originates  from  a solar  dipole  field 
and  not  just  the  effect  of  a magnetic  dipole. 

C.  The  Quadrupole  Terms 

An  examination  of  the  relative  strengths  of  the  dipole  and 
quadrupole  components  of  the  solar  magnetic  field  over  more  than  half 
a solar  cycle  from  1976  through  1983  (Hoeksema,  1984)  reveals  only  two 
instances  in  1980  when  the  magnitude  of  the  quadrupole  component  was 
greater  than  that  of  the  dipole  component,  and  in  each  case  it  was  the 
low  strength  of  the  dipole  which  enabled  the  quadrupole  components  to 
dominate.  Nonetheless,  the  magnitude  of  the  quadrupole  field  was, 


-61- 


except  at  sunspot  minimum,  at  least  30%  of  the  dipole  field  strength, 
and  any  analysis  which  does  not  include  the  effects  of  the  quadrupole 
terms  would  be  incomplete. 

That  part  of  the  double  sum  in  equation  5.4  which  involves  only 
the  quadrupole  terms  can  be  written 

R 4 2 

Bo  = ^ P^coseKg^cosrmj)  + h^sirv^]  . 5.13 

s m= 0 

This  can  be  expanded  out  to  the  five  g and  h coefficients  (i^g  is 
meaningless)  by  substituting  for  P^(cos0)  from  equation  A. 9 to  get 

R 4 i 2 

Bq  = 5(-^)  [^(3cos  e - l)g2g  - 3sinecose(g21cos<|) 

s 

2 

+ h21sin<j>)  + 3sin  ©(g^cosZij)  + i^si^^)]  , 

Rs  2 

and,  since  the  radial  field  drops  off  as  (— ) , 

R 2 R 2 

8 (r)  = 5(^)  (-2)  [i  (3cose  - l)g20 

S 

- 3sinecose(g2jCosi}>  + h^^sin^) 

p 

+ 3sin  e(g22cos2<()  + h22si n24> ) ] . 5.14 

Substituting  for  e and  the  orbital  parameters  relative  to  the  solar 
equatorial  coordinate  system,  i,  fi,  and  u,  we  get 
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R 2 R 2 

Br(r)  = 5(^)  (~)  (3sin2usin2f  - l)g20 

- 3sinusini[(sin^cosu  + cosftsinwcosi)!^ 

+ (cosficosu  - sinfisinucosi^j] 

o 

- 3[(sirU2cosu  + cosfisinucosf) 

2 

- (cosftcosu  - sinftsinucosi)  ]g22 
+ 6(cosftcosu  - sinfisinucosi) 

x (sinficosu  + cosftsinucosi)!^}  • 


The  Lorentz  force  on  the  particle  due  to  the  magnetic  field 
generated  by  a solar  quadrupole  component  is  then 


W = 


5qv  R 2 R 2 , 9 9 

-r~^  (n^)  (-£)  (i  (3sin2usin2f  - 1 


)g 


20 


3sin«sini(si nftcosu  + cosftsinucosi)h 


21 


- 3sinusini(cosftcosu  - sinftsinucosf 

2 

- 3[(sinftcosu  + cosasinucosf) 

2 

- (cosftcosu  - sinfisinucosi ) ]g22 
+ 6(cosftcosu  - sinftsinucosi ) 

x (sinficosu  + cosfisinucosf )h22}  > 


which  affects  only  the  inclination  and  ascending  node  as  was  discussed 
in  the  case  of  the  dipole  field.  These  changes  are 
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u l iiauu  rv  i 


/ 0 \ / 0 \ r 1 / * • 2 ■ 2 

(1  - 3sin  usin  t)g2Q 

''s  ' 


+ 3sinusini(si nficosw  + cosfisi nucosi)!^ 
+ 3sinusini(cosncosu  - sinftsirmcosi^j 


- (cosficosu  - sinfisinucosi)]g22 

- 6(cosftcosu  - sinnsinucosi) 

x (sinftcosu  + cosftsinucosi)!^}  » 


+ 3sinusint(sim2cosu  + cosnsi  nucosi)!^ 
+ 3sinusin£(cosftcosu  - sinftsinucost^j 
+ 3[(sinftcosu  + cosftsi nucosi) 

- (cosftcosu  - sinftsinucosi)]g22 

- 6(cosncosu  - sinftsinucosi) 

x (sinftcosu  + cosnsinucosi)!^}  > 


and 


clt  " nacm 


da  _ Sqv^sinu 


R 2 R 2 

sini  (-^)  (-£)  {|  (1  - 3sin2Ksin2£)g20 


tions 
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5qv  2 2 

^ d i - — 2~^~~  ^g20'^  ^ C0Su  ~ $ sin  ucosu  du] 

(n^acm) 

+ 3cosi[(h21cosfi  - g21sinn)sini  + (g22cos2n 

p 

- h22sin2fi)cosi]  $ sin  ucosu  du 
+3[(h21sinft  + g21cos^)sini  + 2(g22sin2ft 

- h22cos2ft)cosi]  ^ sinucos  u du 

O 

- 3[g22cos2ft  + h22sin2n]  $ cos  u du} 
and 

5qv  3 

<f>  dft  = — 2 — (g?n[.5  ^ sinu  du  - 1.5s ini  $>  sin  u du] 

(n  cicmsini) 

+ 3cosi[(h21cosft  - g21sin^)sini  + (g22cos2n 

O 

- h22sin2^)cosi]  <£  sin  u du 

+3[(h21sinft  + g21cosfi)sini  + 2(g22sin2ft 

2 

- h22cos2fi)cosi]  $ sin  ucosu  du 

2 

- 3[g22cos2ft  + h22sin2ft]  <f>  sinucos  u du} 

are  reached.  On  examining  each  equation,  and  comparing  the  imbedded 
integrals  to  equations  A.1-A.5  in  Appendix  A,  it  can  be  seen  that 

Ai(quadrupole)  = 0 
and 

Afi(quadrupole)  = 0 . 

THE  SURPRISING  RESULT  IS  THAT  THE  QUADRUPOLE  MOMENT  OF  THE  SOLAR 


-65- 


MAGSETIC  FIELD  HAS  NO  EFFECT  ON  THE  ORIENTATION  OF  AN  ORBITING  R AR- 
TICLE'S ORBITAL  PLANE. 


D.  Octupole  Components 

The  octupole  terms  in  the  solar  magnetic  field  can  usually  be 

neglected  since  there  are  only  a few  instances  (in  the  nine  years  of 

observational  record,  1976-1984)  when  the  octupole  field  strength  came 

within  25%  of  the  dipole  or  quadrupole  terms.  Because  a multipole  term 

R 

is  less  by  a factor  of  at  the  source  surface  than  the  next  lower 
order  when  the  two  multipole  orders  are  of  comparable  strength  at  the 
photosphere,  the  higher  order  terms  will  contribute  very  little  to  the 
interplanetary  magnetic  field.  The  octupole  term  is  of  interest 
because  the  lack  of  influence  of  the  quadrupole  field  on  the  orienta- 
tion of  the  particle's  orbit  prevents  the  effects  of  the  octupole 
field  from  becoming  overshadowed. 

The  expanded  expression  for  the  octupole  component  of  the  magnetic 
field  at  the  source  surface  is 

Rq  5 o i 

B0  = 7 (^-)  [P3(cose)g30  + P3(cose) (g31cos<}>  + h31sin<}>) 

p 

+ P3(cose)  (g32cos2(f)  + h32si  n2cj> ) 

3 

+ P3(cose)(g33cos3((>  + h33sin3<j>)  . 

R 2 

We  can  multiply  this  expression  by  (-p)  to  represent  the  field  at 

qv 

the  distance  of  the  particle  and  then  multiply  the  result  by  -p-2  to 
get  the  expression  for  the  force  on  the  particle  due  to  the  octupole 
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component  of  the  magnetic  field.  When  we  then  substitute  the  expres- 
sions for  the  Legendre  polynomials  (equations  A. 9)  we  get 

p 

W = - Fo(r)[g30cose(5cos  e - 3) 

2 2 

- 3g31sino(5cos  e - l)cos<j>  - 3h31sino(5cos  o - l)sin$ 

2 2 2 2 
+ 30g32cosesin  e(cos  <j>  - sin  <j>)  + bOh^cosesin  esin<}>cos<j> 

3 2 2 

- 15g33sin  ecos^cos  $ - 3si n ) 

3 2 2 

- 15h33sin  esi n^> ( 3cos  <}>  - sin  <}>)]  5.15 


where 


F0(D 


3. 5qv 


R 3 R 2 
°'  (-2) 


'R  ' 


It  is  easiest  to  see  the  effects  of  the  individual  component  g's 
and  h's  by  deriving  the  effect  of  each  term  separately  on  i and  n. 

For  g3Q,  after  replacing  those  functions  of  6 and  <t>  with  expressions 
based  on  the  orbital  elements  £,  ft,  and  u,  and  then  substituting  the 
first  term  of  equation  5.15  into  the  Gaussian  perturbation  equations 
(A. 16  and  A. 18),  the  change  in  inclination  after  one  orbit  becomes 


3 3 

Ai(g3o)  = - F^(r)[5sin  i <j>  sin  ucosw  d u 
- 3sin£  f sinucosu  dz^ ] g 3 q 


where 


F_ ( r)  S.Sqv^  R ? R 0 
F,(r)  - 4^  = 4-2  (0,3  (_2j2  ; 
(n  am)  (n  ac)  s 


or 


( 9 sq)  ~ ^ 
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and  the  change  in  the  longitude  of  the  ascending  node  is 

2 4 2 

A^(g30)=  - F^rJtSsin  % <f,  sin  u du  - 3 <f,  sin  u du]g3Q  , 
or 

o 

An(g30)=  3rr  ( r ) [1  - 1.25  sin  i]g3Q  . 

The  g3Q  term  has  no  effect  on  the  inclination,  but  does  precess  the 

2 4 

orbit.  It  should  be  noted  here  that  if  sin  i - the  orbit  will  no 
longer  precess. 

A similar  analysis  can  be  applied  to  the  g31  term.  The  expres- 
sion for  M is 

2 2 2 

Ai(g3i)  = 3F^(r)[5sin  fcosft  ^ sin  ucos  u d u 

2 2 3 

- cosft  § cos  u du  - 5sin  icosfsinft  <£  sin  ucosu  d u 

+ cosisinft  ^sinucosu  du]g3^  . 

When  the  values  for  the  integrals  (equations  A.6-A.8)  are  substituted 
and  the  terms  collected,  the  resulting  expression  becomes 

2 

Ai(g31)  = 3tt  F1(r)cosfi[1.25  sin  i - l]g31  • 

2 4 

When  sin  i = and  cosft  is  negative,  any  increase  in  i which  causes  a 
corresponding  increase  in  the  absolute  value  of  sini  will  cause  a cor- 
responding decrease  in  Ai,  so  that  a stable  resonance  is  reached. 

This  occurs  when  i is  about  63.4  degrees. 

The  expression  for  the  change  in  the  longitude  of  the  ascending 
node  due  to  the  component  of  the  solar  magnetic  field  is 


-68- 


3Fl(r)  2 3 

An(g31 ) = ' si'n^f ' [5sin  icosft  <£  sin  ucosw  d u 

t 2 4 

- cosft  f smucosu  d u - 5s i n fcosfsinft  <J>  sin  u d u 

2 

+ cosfsinft  (f  sin  u dw]g31  , 

which  becomes 

2 

Aft(g31)  = 3 it  F^r)  cotfsinft  [1  - 3.75  sin  f]g31  . 

There  is  also  a resonance  here,  for  when 

cotf[l  - 3.75  sin2 3f]  5.16 

is  positive,  ft  will  proceed  toward  180  degrees,  at  which  point  it  will 
become  stable.  This  independent  analysis  of  the  equilibrium  points 
of  i and  ft  does  not  tell  the  entire  picture,  since  the  stable  point  of 
inclination  results  in  expression  5.16  becoming  negative  which  causes 
ft  to  move  toward  0°.  Then  when  cosft  becomes  positive  (when  ft  reaches 
the  first  or  fourth  quadrants),  the  expression  for  M becomes  positive 
and  any  small  difference  in  i from  the  equilibrium  position  will 
cause  it  to  move  farther  away  from  equilibrium.  We  expect  the  result 
of  this  coupling  of  i and  ft  to  be  a constant  oscillation  between  the 
equilibrium  values  of  i and  ft. 

The  g^  term  in  the  solar  magnetic  field  generates  a change 

Ai ( 9^2 ) = “ ^OF^r)  sini  [(cos2ft  - sin2ft)  $ sinucos  u d u 

2 o 

- (cos2ft  - sin2ft)cos  f ^ sin  ucosu  d u 

2 2 

- 2cosisin2ft  $ sin  wcos  u dujg^  » 
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which  reduces  to 


(932 ) = " *Vr)  sin2fsin2ft  . 

The  corresponding  change  in  the  longitude  of  the  ascending  node  is 

Aft(g32)  = - 30Fj (r)[cos2ft  ^ cos2usin2u  d u 

c • 4 3 

- cos2ft  <j>  sin  u du  - 2cosisin2ft  <f  sin  ucosu  du]g32 

or 

Afi(g22)  = - 7. 5ir  F^(r)  cos2ft(l  - 3cos 2f)  . 

Because  both  expressions  (for  Ai  and  Aft)  have  an  ft  dependence  which 
goes  as  2ft,  any  results  for  a given  ft  will  also  apply  to  an  ft  180 
degrees  greater  (or  smaller).  There  should  be  a stable  equilibrium 
at  2ft  = 270°  and  U = 360°. 

Finally,  the  integral  expressions  for  g33  are 

Ai(g33)  = SOF^r)  {cosft(l  - 4sin2ft)  cos2u  d u 

2 2 2 

- 3cos  f § cos  usin  u d u] 

OOO 

- sinft(l  - 4cos  ft)  [cos  i <f>  sin  ucosk  d u 

- 3cosi  ^ sinucos  u du]}g33  , 

and 

Aft(g33)  = SOF^r)  {cosft(l  - 4sin2ft)  cos\sinu  d u 

2 2 

- 3cos  i § cosusin  u d u] 

2 O A 

- sinft(l  - 4cos  ft)  [cos  i § sin  u d u 

o 9 Q o o 

- 3cosf  & sin  wcos  u du]}  -4^  , 

7 J Qini  5 


-70- 


which  reduce  to 

2 

Af(g33)  = 22 . 5tt  sin  % cos3ft  g^  , 5.17 

and 

( 933 ) = sin2isin3ft  g^  . 5.18 

This  set  of  equations  also  display  equilibrium  positions  for  i and  ft. 

The  inclination  should  reach  a metastable  equilibrium  when  i = 0°  or 

2 

180°  due  to  the  sin  i term  in  the  expression  for  hi.  The  equilibrium 

2 

is  not  stable  since  sin  i is  positive  on  both  sides  of  its  zero  point. 

ft  should  be  stable  when  it  reaches  these  same  values  because  of  the 

sin3ft  factor  in  the  expression  for  Aft. 

The  corresponding  expressions  involving  the  h^  terms  will  not  be 

derived  since  their  only  difference  from  the  g^  terms  is  that  each  is 
90° 

rotated  in  the  <f>  direction  with  respect  to  the  corresponding  g 
term. 


E.  Numerical  Simulations 

In  order  to  simulate  the  orbital  evolution  of  a particle  influ- 
enced by  the  solar  magnetic  field,  the  equations  of  motion  were  solved 
using  Cowell's  Method,  where  the  particle's  position  and  velocity  were 
recomputed  one  thousand  steps  per  orbit.  Table  5.1  shows  the  errors 
in  computation  after  1000  orbits  for  a particle  without  charge, 
initially  in  a circular  orbit  at  3.5  Rq  and  not  undergoing  P-R  drag. 
The  large  error  in  the  argument  of  perihelion  is  meaningless,  since 
there  is  no  perihelion  in  a circular  orbit. 
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Table  5.1.  Numerical  errors  in  the  computational  method. 


el ement 

a (cm) 

e 

i (deg) 

ft  (deg) 

<*>  (deg) 

error 

.06 

4 x 10”7 

6 x 10‘14 

2 x 10'13 

184.5 

It  has  been  shown  by  Lamy  (1976)  that  the  semimajor  axis  remains 
nearly  constant  over  more  than  100  orbits.  This  is  caused  by  a 
balancing  of  P-R  drag  and  the  sublimation  effects  and  results  in  an 
enhancement  in  particle  density  at  the  outer  edge  of  the  dust-free  zone 
surrounding  the  sun.  We  have  used  this  balance  in  order  to  improve 
the  speed  of  computation  by  setting  the  radiation  pressure  coeffi- 
cient equal  to  zero,  which  eliminates  both  the  P-R  drag  and  the  outward 
pseudo-force  mentioned  in  the  introduction.  This  approximation  will 
accurately  portray  the  behavior  of  the  circumsolar  dust  when  the 
interaction  of  the  effects  of  sublimation  and  the  magnetic  field  can 
be  neglected.  This  is  true  to  first  order  since  the  magnetic  field 
does  not  significantly  affect  the  particle-sun  distance,  which 
determines  the  sublimation  rate. 

The  charge  of  the  particle,  q,  admittedly  a subject  of  much 
speculation,  for  purposes  of  this  study  was  set  at  10  volts  (6.7  x 10~^ 
e.s.u.;  cf.  Consolmagno,  1980).  It  is  standard  practice  to  describe 
the  charge  of  a dust  grain  in  volts,  since,  at  least  for  spherical 
particles,  the  grain  potential  is  independent  of  the  size  of  the 
particle.  Since  we  include  particles  of  different  ascending  nodes 
equally  spaced  through  360  degrees,  the  sign  of  the  charge  was 
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determined  to  be  irrelevant  to  the  final  results.  This  is  because  a 
decrease  in  the  (low)  initial  inclinations  of  the  test  particles  causes 
the  reversal  of  the  ascending  and  descending  nodes  and  thus  a particle 
with  a negative  charge  will  experience  the  same  perturbations  as  the 
positively  charged  particle  with  the  opposing  ft.  Also,  the  magnitude 
of  the  charge  can  affect  only  the  timescale  of  the  variations  of  the 
orbital  elements  so  an  error  in  the  estimate  of  the  charge  on  the 
particle  should  not  affect  our  conclusions. 

In  the  first  set  of  simulations,  particles  were  initially  placed 
in  circular  orbits,  at  an  initial  inclination  of  2 degrees  and  given 
six  different  longitudes  of  the  ascending  node  spaced  60  degrees  apart. 
The  orbital  radius  was  set  to  be  on  or  near  the  source  surface,  in 
order  to  simplify  our  understanding  of  the  results,  since  the  simula- 
tion program  was  designed  to  take  into  account  the  finite  speed  of  the 
solar  wind  and  thus  determine  the  strength  of  the  magnetic  field  based 
not  on  the  particle's  angular  coordinates,  but  on  the  position  where 
it  would  have  had  to  leave  the  source  surface  in  order  to  arrive  at 
the  particle's  position. 

Each  such  particle  was  then  given  100  or  more  orbits  in  a mag- 
netic field  consisting  of  every  single  multipole  component  studied 
above.  Thus  for  each  g^  or  h^  , we  have  100  orbits  of  each  of  six 
particles  starting  at  various  values  of  the  longitude  of  the  ascending 
node.  These  results  were  then  compared  to  the  analytical  results.  We 
expect  some  differences  in  the  two  methods  since 

a.  the  analytical  expressions  assume  that  the  non-radial  com- 
ponent of  the  interplanetary  field  in  the  F corona  can  be 
ignored. 


-73- 


b.  the  integration  of  the  expressions  for  M and  am  hold  i and 
M constant  over  one  orbit. 

c.  the  substitution  of  i,  fi,  and  u for  e and  <j>,  which  was  used 
in  the  derivations,  does  not  take  into  account  the  finite 
travel  time  of  the  magnetic  field  from  the  source  surface  to 
the  particle's  position. 

d.  the  linear  step-by-step  nature  of  the  integrations  used  in 
the  simulations  will  introduce  some  error  into  the  computer 
resul ts. 

The  results  of  these  simulations  are  shown  in  figures  5.1  through 
5.10.  The  first  three  figures  are  the  results  for  the  dipole  terms  of 
the  solar  magnetic  field.  Figures  5.4  through  5.6  comprise  the  simu- 
lation results  for  the  g terms  of  the  quadrupole  field,  and  the  last 
four  sets  of  figures  show  the  variations  in  i and  m for  the  g terms  of 
the  octupole  field.  These  computer  simulations  were  especially  useful 
in  studying  the  equilibrium  points  of  i and  m generated  by  the  octupole 
field  terms.  The  corresponding  h terms  have  been  suppressed  since 
there  exists  only  a rotation  equivalence  between  each  g term  and  its 
corresponding  h term. 

Figure  5.1a  demonstrates  that  the  inclination  did  not  change  for 
any  of  the  six  (superimposed)  starting  values  of  the  longitude  of  the 
ascending  node  when  a polar  dipole  field  was  modelled.  The  linear 
form  of  the  variations  of  the  longitudes  of  the  ascending  nodes  vs. 
time  in  figure  5.1b  corresponds  to  the  relation  in  equation  5.12b. 
Rather  than  calculate  the  rate  of  change  of  the  longitude  of  the 
ascending  node  from  equation  5.12b,  we  can  determine  it  directly  from 
figure  5.1b,  as  -1.1  degrees  per  orbit. 
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(a) 


60  80  100  120  140  160  l8«  200 

NUHBER  OF  ORB  I T3 


(b) 


Figure  5.1.  The  effect  of  a g-iQ  (polar)  dipole  field  on  a)  inclina- 
tions and  b)  the  longitudes  of  the  ascending  nodes  of 
dust  orbiting  near  the  sun. 
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The  inclination  change  due  to  the  component  of  the  solar 
magnetic  field  is  constant  with  respect  to  both  the  inclination  and 
the  longitude  of  the  ascending  node,  as  seen  in  figure  5.2a.  Observe 
that  the  ascending  nodes  of  the  six  particles  in  figure  5.2b  all  quickly 
align  themselves  to  fi  = 0°,  until  the  inclination  reaches  180°  after 
approximately  170  orbits,  at  which  point  the  ascending  and  descending 
nodes  flip,  and  the  longitude  of  the  ascending  node  becomes  180°. 

The  rate  of  change  of  inclination  is  then  slightly  less  than  1°  per  orbit. 

As  has  been  stated  above,  the  hll  components  of  the  solar  magnetic 
field  have  the  same  form  as  the  corresponding  g^  components,  and  the 
only  difference  in  their  influence  is  a rotation  in  the  e direction; 
hence,  there  will  be  a difference  in  any  alignment  of  the  ascending 
nodes.  Figures  5.3a  and  5.3b  demonstrate  this  fact.  Here  the  effect 
on  inclination  is  the  same  as  that  for  figure  5.2a,  but  the  curves  of  U 
in  figure  5.3b  are  displaced  90°  from  those  in  figure  5.2b. 

Figures  5.4,  5.5,  and  5.6  come  as  no  surprise  since  we  have 
already  shown  that  the  quadrupole  components  have  no  effect  on  the 
orientation  of  the  particle's  orbit.  It  is  gratifying  to  have  a null 
result  in  the  simulations,  but  more  importantly,  these  results  verify 
that  the  approximations  used  in  deriving  the  analytical  expressions 
are  valid. 

Once  again  the  difference  between  the  simulations  and  the  analy- 
tical derivations  must  be  emphasized.  The  computer  simulations  were 
performed  by  calculating  the  acceleration  due  to  the  magnetic  field 
(and  the  sun)  at  a time  interval  of  one  thousandth  of  an  "orbital 
period,"  and,  by  using  these  accelerations, computing  thenew  position 
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Figure  5.2.  The  effect  of  a g-Q  dipole  field  on  a)  inclinations  and 

b)  the  longitudes  of  the  ascending  nodes  of  dust  orbiting 
hear  the  sun. 
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Figure  5.3.  The  effect  of  an  hn  dipole  field  on  a)  inclinations 
and  b)  the  longitudes  of  the  ascending  nodes  of  dust 
orbiting  near  the  sun. 
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(a) 


(b) 


Figure  5.4.  The  effect  of  a g20  quadrupole  field  on  a)  inclinations 
and  b)  the  longitudes  of  the  ascending  nodes  of  dust 
orbiting  near  the  sun. 
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(b) 


Figure  5.5.  The  effect  of  a g2i  quadrupole  field  on  a)  inclinations 
and  b)  the  longitudes  of  the  ascending  nodes  of  dust 
orbiting  near  the  sun. 
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Figure  5.6.  The  effect  of  a g2?  quadrupole  field  on  a)  inclinations 
and  b)  the  longitudes  of  the  ascending  nodes  of  dust 
orbiting  near  the  sun. 
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and  velocity  vectors  of  the  particle  at  that  time.  This  method 
neither  uses  the  particle's  orbital  elements  (except  to  redetermine 
the  time  interval  between  subsequent  calculations  after  every  print- 
out) nor  does  it  integrate  the  effects  over  one  orbit.  Indeed,  the 
computations  do  not  presuppose  the  existence  of  an  osculating  orbit, 
although  such  orbital  elements  are  generated  for  the  output,  since 
x,  y,  z values  are  not  useful  for  diagnostic  purposes. 

The  results  for  the  effect  of  the  octupole  field  components  on 
the  orientation  of  orbiting  dust  particles  are  by  far  the  most  exciting. 
Because  the  octupole  components  are  relatively  weak  throughout  most 
of  a solar  cycle,  it  is  not  expected  that  these  results  will  influence 
the  structure  of  the  F coronal  dust  cloud  to  a great  extent.  However, 
due  to  the  fact  that  most  of  the  expressions  derived  in  section  D 
indicate  that  each  octupole  term  creates  stable  or  metastable  orbital 
alignments,  it  is  possible  that,  at  times  of  a strong  octupole  term, 
there  would  be  a stable  structure  created  temporarily  in  this  region 
which  could  be  observed.  This  suggests  conversely  that  such  an 
observation  would  confirm  the  presence  of  a relatively  strong  octupole 
field,  which  should  be  of  interest  to  solar  physicists. 

Figure  5.7a  demonstrates  the  lack  of  effect  of  the  g^Q  component 
on  the  inclinations  of  orbiting  particles.  The  straight  line  across 
the  bottom  of  the  graph  represents  the  inclinations  for  all  six  of  the 
initial  values  of  the  longitude  of  the  ascending  node.  Since  ft  varied 
with  time  (figure  5.7b)  in  all  six  cases,  while  i did  not,  there  appears 
to  be  no  relation  between  M and  ft.  In  fact,  this  result  is  identical 
with  that  in  section  D.  From  figure  5.7b,  we  can  determine  that  the 
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Figure  5.7 


The  effect  of  a g3Q  octupole  field  on  a)  inclinations 
and  b)  the  longitudes  of  the  ascending  nodes  of  dust 
orbiting  near  the  sun. 
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change  in  the  longitude  of  the  ascending  node  is  about  120°  for  200 
orbits  or  0.6°  per  orbit. 

As  described  in  section  D,  the  g^  octupole  component  drives  the 
inclination  toward  about  63°,  which  can  be  seen  in  figure  5.8a.  The 
longitude  of  the  ascending  node  quickly  reaches  an  equilibrium  value 
of  180°  in  the  case  of  all  six  particles  starting  at  differing  values 
for  ft.  The  reason  for  the  rapid  change  in  ft  compared  to  i is  due  to 
the  coti  term  in  the  expression  for  a ft,  since  at  the  beginning  of 
these  simulations  the  inclination  is  2°,  which  makes  coti  very  large. 
Notice  also  how  ft  diverges  from  180°  when  i becomes  greater  than  about 
31°,  the  point  where  sin^i  becomes  y^. 

The  stable  equilibrium  value  for  ft  due  to  the  g32  term  is  at  2ft 
= 90°  or  ft  = 45°  or  225°,  as  can  be  seen  in  figure  5.9b.  This  causes  i 
to  proceed  toward  90°,  where  sin2i  becomes  0 (figure  5.9a).  Once  cos 2i 
becomes  less  than  y,  the  expression  for  the  change  in  the  longitude  of 
the  ascending  node  becomes  negative,  and  the  value  of  ft  diverges  from 
the  equilibrium  positions. 

In  figures  5.10a  and  5.10b,  we  can  see  the  effect  of  the  g33  term 
on  the  orientation  of  the  particle's  orbit.  As  expected,  the  sin3ft 
term  in  the  expression  for  Aft  (equation  5.18)  causes  the  change  in  ft  to 
become  small  when  ft  is  near  an  integral  multiple  of  60  degrees.  The 
reason  the  values  which  give  3ft  = 180°  have  larger  variations  is  that 
the  sign  of  sin3ft  drives  3ft  away  from  180°  and  toward  0°  or  360°.  The 
variation  in  inclination  also  displays  the  threefold  symmetry  in  the 
longitude  of  the  ascending  node,  so  that  only  two  curves  of  i vs.  time 
are  visible.  The  inclination  change  is  small  due  to  the  sini  term  in 
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Figure  5.8 
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(a) 


Figure  5.9.  The  effect  of  a 930  octupole  field  on  a)  inclinations 
and  b)  the  longitudes  of  the  ascending  nodes  of  dust 
orbiting  near  the  sun. 
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Figure  5.10.  The  effect  of  a 933  octupole  field  on  a)  inclinations 
and  b)  the  longitudes  of  the  ascending  nodes  of  dust 
orbiting  near  the  sun. 
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the  numerator  of  the  expression  for  i (equation  5.17).  Notice  also 
that  the  change  in  i accelerates  when  i becomes  larger,  due  to  the 
same  sini  term. 

A second  type  of  simulation  was  run  based  on  the  observed  solar 
magnetic  field.  We  believe  it  Is  important  to  determine  the  effect 
of  the  solar  magnetic  field  when  the  entire  field  structure  is  taken 
into  account.  This  type  of  modelling  allows  any  interactions  between 
the  multipole  components  to  be  observed,  as  well  as  determining  what 
sort  of  structure  in  the  solar  dust  cloud  can  be  generated  in  the  F 
coronal  region.  These  results  can  (and  will)  be  compared  with  those 
infrared  observations  of  the  F corona  which  were  briefly  mentioned  in 
the  introduction. 

The  model  of  the  solar  magnetic  field  used  in  this  section  is 
based  on  actual  observations  of  the  solar  surface  and  the  inter- 
planetary magnetic  field.  We  have  used  the  numerical  representation 
of  the  field  developed  by  Hoeksema  (1984)  in  which  the  field  is 
generated  from  the  various  multi  pole  components  of  coefficients  g^m 

and  h.  . These  coefficients  were  determined  from  observations,  with 

Lm 

the  assumption  of  a spherical  source  surface.  In  our  calculation  of 
the  solar  field,  we  used  only  terms  up  to  and  including  the  second 
order  (quadruopole)  coefficients.  Since  our  particle  orbits  are  be- 
yond the  source  surface,  determination  of  the  field  at  any  position 
involves  calculating  the  field  at  that  position  on  the  source  surface 
where  it  originated  and  projecting  it  outward,  taking  into  account 
the  finite  velocity  of  the  field  and  the  effects  of  solar  rotation. 
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Particles  were  run  at  initial  inclinations  of  2°,  5°,  10°,  and 
30°,  and  at  20  different  ascending  nodes  from  18°  to  360°  inclusive. 
Each  case  was  run  at  various  starting  times  so  that,  at  the  date  under 
study,  particles  would  have  been  in  the  region  from  0 to  90  orbits. 
Each  particle  thus  was  represented  by  a set  of  orbital  elements.  The 
x,  y,  z coordinates  of  up  to  twenty  points  (based  on  the  fan  model  of 
Leinert,  1976,  although  our  results  are  not  model  dependent)  along 
each  orbit  were  computed.  These  points  were  equally  spaced  in  the 
mean  anomaly  with  a random  starting  position  of  the  first  point  in 
order  to  avoid  particle  alignments  due  to  the  choice  of  positions. 

The  points  were  then  projected  onto  a viewing  plane  and  plotted.  Fig- 
ure 5.11  shows  the  initial  particle  positions,  including  all  inclina- 
tions and  ascending  nodes,  for  one  starting  time.  We  made  three  runs 
at  different  dates  to  represent  different  configurations  of  the  solar 
magnetic  field.  1976.6  represents  the  field  at  sunspot  minimum,  where 
the  polar  dipole  component  strongly  dominates  the  higher  orders, 

1980.6  (sunspot  maximum)  is  a time  when  the  dipole  field  is  weak  and 
the  quadrupole  terms  dominate,  and,  around  1982.9,  the  contributions 
from  the  dipole  and  quadrupole  components  are  comparable.  Figures 
5.12,  5.13,  and  5.14  show  the  direction  and  magnitude  of  the  field  in 
solar  equatorial  coordinates,  for  1976.6,  1980.6,  and  1982.9,  respec- 
tively. Figures  5.15,  5.16,  and  5.17  show  the  equatorial  dipole 
components  of  these  fields. 

Figures  5.18  and  5.19  show  the  structure  of  the  solar  dust  cloud 
during  sunspot  minimum.  The  particles  have  spread  from  their  initial 
concentrations  (figure  5.11)  to  create  a nearly  homogeneous  band  of 
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Figure  5.11.  Input  particle  distribution  for  the  simulations  of  F 
coronal  dust  distribution  under  the  effect  of  a solar 
magnetic  field. 
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Figure  5.18.  Particle  distribution  in  the  F corona  for  the  date 
1976.6,  as  viewed  from  the  -y  direction.  The  x 
direction  is  the  direction  of  0°  solar  longitude  at  the 
beginning  of  a Carrington  Rotation.  The  distribution 
combines  particles  released  at  times  from  0 to  90  orbits 
before  1976.6 
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Figure  5.19.  Particle  distribution  in  the  F corona  for  the  date 
1976.6,  as  viewed  from  the  +x  direction. 
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material  up  to  25  degrees  above  and  below  its  symmetry  plane,  which 
is  no  longer  aligned  with  the  solar  equatorial  plane,  with  particle 
orbits  up  to  40°  from  the  symmetry  plane.  Since  the  solar  magnetic 
field  at  sunspot  minimum  approximates  a polar  dipole,  this  field  was 
not  expected  to  cause  the  radical  increases  in  inclination  displayed 
at  other  times. 

In  figures  5.20  and  5.21,  the  solar  dust  zone  has  been  under  the 
influence  of  the  quadrupole  dominated  solar  magnetic  field  of  1980.6. 
Once  again,  the  particles  have  been  moved  from  their  original  orbits 
so  that  the  plane  of  symmetry  is  no  longer  the  solar  equatorial  plane. 
The  structure  in  figure  5.11  is  evidence  of  an  aligning  of  the  lines 
of  nodes  of  the  particle  orbits  preferentially  in  the  direction  of  the 
x-axis.  Comparing  this  with  the  dipole  component  of  the  field  in 
1980.6  (figure  5.16),  it  appears,  as  expected,  that  the  orbital  planes 
have  precessed  about  the  line  of  the  maximum  dipole  component. 

Figures  5.22  and  5.23  demonstrate  that  the  solar  magnetic  field 
can,  at  times,  completely  scatter  the  particles  from  the  solar  equa- 
torial plane.  The  inclinations  of  the  particle  orbits  go  from  0°  to 
120°,  where  inclinations  of  greater  than  90°  apply  to  particles  in 
retrograde  orbits  (Shcheglov  et  al.,  1985,  observed  a retrograde  com- 
ponent in  their  radial  velocity  measurements).  There  is  no  evidence 
of  ring  or  band  structure;  the  distribution  of  the  solar  dust  cloud 
appears  as  a spherical  shell  from  the  direction  of  the  x-axis  and 
covers  about  of  this  region  when  viewed  along  the  y-axis.  Once 
again  the  equatorial  dipole  field  is  the  dominant  factor  in  the  spread 
of  particle  inclinations,  as  can  be  seen  by  comparing  the  axis  of 
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Figure  5.20.  Particle  distribution  in  the  F corona  for  the  date 
1980.6,  as  viewed  from  the  -y  direction. 
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Figure  5.21.  Particle  distribution  in  the  F corona  for  the  date 
1980.6,  as  viewed  from  the  +x  direction. 
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Figure  5.22.  Particle  distribution  in  the  F corona  for  the  date 
1982.9,  as  viewed  from  the  -y  direction. 
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Figure  5.23.  Particle  distribution  in  the  F corona  for  the  date 
1982.9,  as  viewed  from  the  +x  direction. 


-103- 


intersection  of  the  orbital  planes  with  the  equatorial  dipole  field 
in  figure  5.17.  The  results  of  Hoeksema  (1984)  indicate  that,  for  the 
period  1976  to  1984,  the  equatorial  dipole  component  was  as  strong  or 
stronger  than  that  of  1982.9  about  half  of  the  time. 

Figure  5.24  and  5.25  show  the  distribution  of  just  those  particles 
with  initial  inclinations  of  2°.  The  band  structure  of  each  group  of 
particles  with  the  same  starting  time  remains  nearly  intact,  but  each 
such  band  has  been  reoriented  in  space.  A study  done  on  the  change 
in  the  particles'  orbital  planes  to  resolve  the  variation  in  the  pre- 
cession axis  determined  that  the  axis  did  not  change  more  than  5°  in 
orientation  over  the  course  of  ninety  orbits. 

F.  Analysis  of  Observational  Evidence 

The  observations  mentioned  in  the  introduction,  which  have  shown 
an  enhancement  at  approximately  4RQ,  are  all  within  30°  of  the  ecliptic 
plane.  Unfortunately,  those  observations  attempted  by  MacQueen  et 
al.  that  would  have  scanned  perpendicular  to  the  ecliptic  were  lost 
due  to  equipment  failure,  so  that  there  is  no  solid  evidence  as  to  the 
maximum  extent  of  the  enhancement,  although  there  appears  to  be  no 
limit  in  the  latitude  at  which  the  dust  appears.  Previous  theoreti- 
cians have  described  the  enhancement  in  density  as  a ring  (Mukai , 

1983).  They  based  their  belief  on  the  observation  that  the  inter- 
planetary zodiacal  cloud  indicates  a much  higher  dust  density  within 
10°of  the  invariable  plane.  However,  a dust  ring  aligned  with  any 
plane  not  normal  to  the  ecliptic  will  appear  from  the  Earth  as  an 
ellipse,  a thin  ellipse  if  the  viewer  is  within  a few  degrees  of  the 
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Figure  5.24.  Particle  distribution  in  the  F corona  for  the  date 
1982.9,  as  viewed  from  the  -y  direction.  Only 
particles  with  initial  inclinations  of  2°  are 
included. 
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Figure  5.25.  Particle  distribution  in  the  F corona  for  the  date 
1982.9,  as  viewed  from  the  +x  direction.  Only 
particles  with  initial  inclinations  of  2°  are 
included. 
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ring  plane,  for  the  maximum  length  of  the  minor  axis  of  such  an 
ellipse  as  viewed  from  Earth  is  .5  RQ  when  the  Earth  is  90°  or  270° 
from  the  direction  of  the  line  of  nodes  of  the  solar  equatorial  plane. 

Table  5.2  shows  a summary  of  the  observations.  Here,  XQ-ft0  is 
the  angle  between  the  Earth  and  the  ascending  node  of  the  solar  equa- 
torial plane  at  the  time  of  observation.  Scan  angle  I is  the  angle 
between  the  line  of  observation  and  the  ecliptic,  and  scan  angle  II 
(calculated  by  the  author)  is  the  angle  between  the  line  of  observation 
and  the  solar  equatorial  plane,  assuming  spherical  symmetry.  While 
there  are  no  observations  at  high  solar  latitudes,  there  is  enough 
latitude  variation  in  the  available  observations  to  show  that  there 
cannot  be  a structure  in  the  solar  dust  cloud  which  does  not  cover  22° 
north  or  south  of  the  solar  equatorial  plane,  and  this  is  a minimum. 

Numerical  simulations  of  the  effect  of  the  solar  magnetic  field 
on  the  distribution  of  the  coronal  dust  particles  show  that  the  mag- 
netic field  tends  to  scatter  particles  toward  the  polar  regions,  with 
a significant  effect  arising  when  the  equatorial  dipole  components 
become  strong.  According  to  the  observational  evidence  (Hoeksema, 

1984,  p.  177),  this  is  at  least  50%  of  the  time.  Observations  appear 
to  support  the  conclusion  that  the  dust  cloud  exists  as  a wide  band  or 
spherical  shell  surrounding  the  sun.  Further  observational  studies 
of  infrared  emission  in  the  F corona,  especially  at  high  solar  lati- 
tudes, are  needed  to  verify  these  conclusions. 
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Table  5.2.  A summary  of  observations  of  the  enhancement  in  the  infrared 
emission  at  4 R0.  Scan  angles  I and  II  are  measured  from 
the  line  of  observation  to  the  ecliptic  plane  and  the  solar 
equatorial  plane,  respectively.  Scan  angles  II  were  cal- 
culated by  this  author. 
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Date 
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Angle 
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Observed 

Feature 

r/R0 

Scan 

Angle 

II 

Beavers  et  al . 

6/26/79 

98° 

n/a 

4.3 

21° 

MacQueen 

11/12/66 

203° 

0° 

3.9 

7° 

MacQueen 

1/  9/67 

146° 

20° 

4.2 

14° 

MacQueen 

1/  9/67 

146° 

-40° 

4.2 

46° 

Maihara  et  al . 

6/11/83 

355° 

7° 

4. 

14° 

Peterson 

9/  3/66 

272° 

? 

M. 

- 

Peterson 

11/12/66 

203° 

29° 

4.2 

22° 

Peterson 

3/  7/70 

89° 

var. 

4. 
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CHAPTER  VI 


THE  EFFECT  OF  THE  SUN'S  FINITE  SIZE  AND  SOLAR  ROTATION 


In  the  region  of  the  F corona,  solar  radiation  reaches  an  orbiting 
particle  from  a wide  range  of  directions  due  to  the  fact  that  the  sun 
is  not  a point  source.  Indeed,  a point  just  above  the  solar  photosphere 
receives  radiation  as  if  the  source  were  an  infinite  plane.  Because 
both  the  differential  source  surfaces  and  the  orbiting  particle  are 
moving  with  respect  to  the  center  of  the  sun,  radiation  from  different 
directions  will  impinge  upon  the  particle  with  different  average 
energies.  An  anisotropy  in  the  radiation  field  is  thus  created, 
which  one  must  take  into  account  when  the  effect  of  solar  radiation 
on  the  dynamics  of  an  orbiting  particle  is  to  be  determined. 

In  1962,  Guess  derived  an  expression  for  the  force  generated  by 
radiation  from  a non-rotating  spherical  source  using  the  physics  of 
special  relativity.  His  derivation  did  not  utilize  Robertson's  1937 
expression  for  the  radiation-induced  force  from  a radiating  point 
source,  and  in  fact,  Guess  followed  the  lines  of  Robertson's  deriva- 
tion, essentially  rederiving  the  P-R  drag  formula  in  the  process.  As 
a by-product  of  our  derivation,  we  will  show  that  his  expression  can 
be  derived  without  using  a relativistic  approach.  Here  I also  wish 
to  note  that  although  Robertson's  formula  was  determined  from  rela- 
tivistic considerations,  it  too  can  be  derived  classically  (see  Burns, 
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Lamy,  and  Soter,  1979).  In  fact,  all  of  the  radiation  pressure  effects 
on  orbiting  bodies  can  be  treated  as  classical  forces,  since  the 
velocities  in  question  are  much  much  smaller  than  the  speed  of  light. 

The  forces  on  a particle  due  to  the  finite  size  of  the  sun  and 
due  to  solar  rotation  can  be  derived  jointly,  since  both  involve  an 
integration  of  the  incremental  forces  of  radiation  across  the  solar 
surface.  In  this  section  we  intend  to  determine  the  radiation  force 
from  each  source  point  by  determining  the  relative  velocity  of  the 
source  and  particle  and  apply  Robertson's  (1937)  expression  for  the 
force  due  to  a point  source  sun  at  each  such  source  point. 

We  will  set  up  our  primary  coordinate  system  such  that  the  particle 
is  on  the  x-axis  at  a distance  r from  the  sun's  center,  the  y-axis  is 
perpendicular  to  the  solar  rotation  axis  and  positive  in  the  direction 
of  motion  of  the  particle's  subsolar  point,  and  the  z-axis  is  such 
that  x-y-z  form  the  usual  right-handed  coordinate  system. 

Because  the  velocity  of  a point  on  the  solar  surface  depends  on 
its  position  with  respect  to  the  solar  equatorial  plane,  we  will  need 
to  define  a solar  coordinate  system  xQ,  y , zQ,  such  that  the  zQ-axis 
is  in  the  direction  of  the  sun's  rotation  vector  and  the  yQ-axis  is 
parallel  to  the  y-axis.  The  transformation  between  the  x-y-z  coordinates 
and  the  x0-y0-z0  coordinates  can  be  written  as  a function  of  the  angle 
9,  the  angle  between  the  solar  north  pole  and  the  x-axis.  They  are 

xQ  = x sine  - z cose 

y0  - y 

zQ  = x cose  + z sine  . 
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If  we  consider  a point  S on  the  solar  surface,  with  solar  coor- 
dinates xQ,  y , zQ,  the  velocity  vector  of  that  point  is 


The  velocity  of  the  particle  in  this  system  does  not  depend  on  the 
source  position  and  will  be  denoted 


The  relative  velocity  of  the  particle  with  respect  to  the  source  point 
is 


For  the  purpose  of  determining  the  contribution  of  each  source 
area,  it  is  simplest  to  make  another  transformation  to  a circular 
disk,  i.e. , the  apparent  solar  disk  as  viewed  from  the  particle.  This 
is  done  by  projecting  each  source  point  S onto  the  disk  (see  figure 
6.1)  at  a point  S'.  The  point  S'  has  coordinates  denoted  by  X,Y,Z. 

From  the  geometry  of  figure  6.1,  we  can  determine  the  radius  of  the  disk 
RD,  and  X,  which  is  the  distance  from  the  disk  to  the  center  of  the 
sun.  The  distance  from  the  particle  to  the  edge  of  the  sun  along  a 
tangent  (distance  PQ  in  figure  6.1)  is 


s 


or,  in  our  base  coordinates. 


v$  = -to0y  sine  x + u>Q(x  sine  - z cose)  y + ioQy  cose  z . 


vp  = (x,y,z)  = x x + y y + z z . 
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Figure  6.1.  The  geometry  of  radiation  from  the  solar  surface.  P is 
the  particle  position,  the  source  of  the  surface  is  S, 
and  the  projection  of  S onto  the  apparent  solar  disk  is  S'. 


and  so  the  radius  of  the  disk  is  given  by 


rd2  = r2  - RQ2  - (r  - X)2 


which  also  equals 


Rq2  = R02  - X2 


R 2 R 2 

This  gives  X = and  so  R^2  = Rq2 ( 1 - 


). 


i 

Now,  the  relation  between  the  coordinates  x,y,z  of  S and  the 
coordinates  X,Y,Z  of  S'  is  such  that,  if  we  define  an  angle  e between 
the  line  PS  and  the  x-axis, 


tanE.£HZ.£LlZ 

(r  - x)  (r  - X) 


or 
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9 


and,  because  of  rotational  symmetry, 


V _ Z _ (r  - X) 
y " z (r  - x)  • 


6.1 


as 


2 2 2 

Using  the  fact  that  x + y + z 


we  can  write  equation  6.1 


r 


(r  - X)y 


6.2 


Since  r and  X are  constants,  we  will  define  r - X e p.  Substituting 
2 = Xj  into  equation  6.2,  we  get 


By  squaring  both  sides,  we  can  solve  for  ^ using  the  quadratic 
formula.  We  get 


[2rP  - \l 4r2p2  - 4rp(p2  + R2)] 
2(p2  + R2) 


6.3 


2 2 2 

where  R =X  +Y.  R is  in  fact  the  radial  distance  from  the  center 
of  the  disk  to  the  point  S' , and  it  will  be  convenient  in  later  inte- 
grations to  use  the  transformation  Y = R cos4> , and  Z = R sin<K  which 

in  no  way  invalidates  our  previous  definition  of  R. 

2 2 2 

After  defining  p + R = D , which  equals  the  distance  from  the 
source  point  to  the  particle,  and  collecting  terms,  equation  6.3  becomes 
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y = 


Y[rP  - JrP(RD2 


D 


2 


We  will  cal  1 


[rp  - Jr p(Rd2  - R2)] 


D 


E f 


Remember  that  f and  D are  functions  of  R.  The  expressions  for  x,  y, 
and  z in  terms  of  X,  Y,  and  Z are 

x = r - fp 
y = fY 
z = fZ  . 


The  coordinates  of  the  source  point  in  our  disk  coordinate  system 
are  (X,Y,Z),  and  the  coordinates  of  the  particle  are  (r,0,0).  The 
radius  vector  from  the  source  to  the  particle  is 

R 2 

r = (r  - — ,-Y ,-Z)  = (p  ,-Y ,-Z) 
r 

and  the  relative  velocity  vector,  in  terms  of  the  disk  coordinates,  is 

vrei  = (x  + (^Yfsine,  y - u)Qrsin0  + u0pfsine  + u^Zfcose, 
z - coQYfcose)  . 

The  radial  velocity,  r,  can  be  found  from  r = (-p)-vre^ 
r = -jj  (px  - Yy  - Zz  + co0rYsin6)  . 


, which  gives 
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Now,  with  r,  r,  and  v in  terms  of  the  disk  coordinates  of  the 
source  point,  we  can  substitute  these  values  into  Robertson's  (1937) 
formula 


using  Eq  here  as  the  momentum  flux  per  unit  surface  area  of  the  visible 
solar  disk  impinging  on  a particle  of  unit  size.  We  then  get 


T ~ x HI  " 

~ ^4  (px  “ 

Yy  - 

Zz  + u0rYsine)  - 

x 1 
2J 

to  DJ 

CD* 

cD^ 

+ 

<<  > 
i — i 

+ 4 (px 

- Yy 

- Zz  + w rYsine) 

- -ij] 

cD4 

0 

cD2 

+ (px 

cD4 

- Yy 

- Zz  + u^rYsine) 

cD2 

By  holding  ui  constant,  this  expression  can  be  integrated  over  the 
visible  solar  disk.  does,  in  fact,  vary  with  latitude,  but  the 
main  contribution  to  the  solar  rotation  induced  perturbations  of  a 
dust  particle's  orbit  is  from  the  equatorial  regions  of  the  sun,  where 
the  rotation  rate  is  nearly  constant.  If  we  substitute  R and  <t>  for  Y 
and  Z,  we  can  eliminate  those  terms  which  include  a functional  depen- 
dence on  <(>  of  sin<j>,  cos<|>,  or  sin^costj),  since  the  integral  (from  0 to 
2tt)  of  these  functions  is  zero.  This  leaves 


D 


~ o w r 

+ y [(\)cos2<}>  [-J-  sine 
D4  c 

- i(|)  [(4)sin2*  ■ 

c D4 


c' 


cD^ 
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The  increment  of  surface  area  on  the  disk  as  "seen"  by  the  particle  is 
(^jy)d<j>dR,  where  ^ is  the  cosine  of  the  angle  between  the  particle- source 
line  and  the  normal  to  the  disk.  Integrating  then  over  the  entire 
disk,  we  get,  for  the  directions  x,  y,  and  z, 


R. 


Fx=*Eo^2-  c 


2x  r4 


R023/2 

2 } 
r 


- (1 


^)1/2]}  . 
r 


6.4a 


» E0  - 3<‘  - -3-> 


V 1/2 

2 


1 Rn2  3/2  ? Rn2 

+ 3 ( 1 - -y-)  ] + [7  - (1  - ~2~) 


1 Rn  3/2 

+ 3 (1  ’ 


6.4b 


F = 


7 ft  Rn  1/2 

Eo  ^ [f  - 3<‘  - 


*3  I1  ' 


R/  3/2 

-jr)  ]} 


6.4c 


respectively.  These  expressions  can  best  be  understood  by  expanding 
R 

them  in  terms  of  These  then  become 
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Fx  =’Ro2 *  -§  U 
r 


2X  r,  . 1 R0 

“ L1  + 24 


+ 1 Ro\ 
+ ^~7  + 


. ]>  . 


6. 5a 


2 R H 

F - d 2 r -y  r,  1 Ro  ,5  — 

Fy  _irRo  ^2  { c [1  + 2 ~^2  + 2*  r4 * 

, R 6 . u R Rn  , R* 

+ +?  (-F><->i™  [i 

o R 4 , R 6 

+ A 9 + A 

16  r4  8 r6J}  ’ 


6.5b 


E * R 2 n 4 

F - n 2 Lo  ,-z  n , 1 Ko  ,5  Ko 

Fz  -Ro  "2  {_^  [1  + 2~2  + 24~4 
r r r 


R 


+ ^“T]}  * 


6.5c 


2 

The  expression  it  Rq  in  equations  6.5a,  b,  and  c is  the  area  of 

2 

the  visible  solar  disk  when  viewed  from  a large  distance,  so  i c RQ  Eq 
is  the  total  energy  emitted  from  a solar  hemisphere  at  unit  distance 
impinging  on  a particle  of  unit  size.  If  the  solar  rotation  is  neg- 
lected, the  expressions  in  the  y and  z directions  become  equivalent 
and  equal  to  the  result  of  Guess  (1962).  The  ratio  of  the  force  due 
to  solar  rotation  to  that  of  P-R  drag  at  large  distances  from  the  sun 

i s 


Fsr  . 1 ,“oRo2, 

Fpr~*(  r ' 


sine 
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If  we  enter  the  values  for  u)Q  and  R , and  assume  the  particle  is  orbit- 
ing in  the  solar  equatorial  plane,  so  that  sine  = 1 and  z = z = 0,  we 
can  derive  the  ratio  at  several  different  distances  from  the  sun. 

This  result  has  been  tabulated  in  table  6.1. 


Table  6.1.  The  ratio  of  the  force  due  to  solar  rotation  compared  to 
the  P-R  drag  for  many  values  of  the  particle's  orbital 
radius. 


Radius 

Fsr 

Fpr 

Radius 

Fsr 

Fpr 

10.0000  AU 

0.0000257 

0.1000  AU 

0.0002574 

8.0000  AU 

0.0000288 

16.0000  RQ 

0.0002982 

6.0000  AU 

0.0000332 

12.0000  RQ 

0.0003442 

4.0000  AU 

0.0000407 

10.0000  RQ 

0.0003769 

2.0000  AU 

0.0000576 

8.0000  RQ 

0.0004210 

1.0000  AU 

0.0000814 

6.0000  R 

0 

0.0004851 

0.8000  AU 

0.0000910 

4.0000  R„ 
0 

0.0005908 

0.6000  AU 

0.0001051 

3.5000  R 

0 

0.0006297 

0.4000  AU 

0.0001287 

3.0000  Rn 

0.0006771 

0.2000  AU 

0.0001821 

1.0000  Rn 
0 

0.0011294 

Notice  that  the  ratio  decreases  as  the  distance  from  the  sun 
increases  and  is  never  large,  even  very  close  to  the  sun.  In  the  F 
corona,  although  the  ratio  is  small,  the  P-R  force  is  balanced  by  the 
sublimation  effect  and  the  force  due  to  solar  rotation  is  not 
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overpowered.  This  is  true  only  in  the  direction  perpendicular  to  the 
orbital  plane  of  the  particle,  since  the  coplanar  component  of  the 
force  due  to  solar  rotation  also  becomes  part  of  the  P-R,  sublimation 
balance.  The  perpendicular  component  of  the  solar  rotation  force 
affects  the  inclination  and  the  longitude  of  the  ascending  node  of 
orbiting  particles.  We  can  find  the  magnitude  of  the  acceleration  from 
this  component  by  taking  the  dot  product  of  the  force  vector  due  to 
solar  rotation  with  the  unit  vector  in  the  direction  of  the  particle's 
orbital  angular  momentum  and  dividing  by  the  mass  of  the  particle,  mg. 
In  the  solar  equatorial  coordinate  system,  and  using  the  approximation 
that  ~ is  small  at  the  evaporation  zone),  this  gives 


or,  when  we  assume  a circular  orbit  where  e = 0,  r = a,  and  d u = 


Substituting  the  expression  for  Fsr  (equation  6.7)  into  equation 


-F 


sinfcosu 


W. 


sr 


6.6 


sr 


m 


9 


where 


6.7 


d 7* 

From  Appendix  A,  equation  A. 16,  the  expression  for  is 


o 

(-=p^)  dt,  we  can  substitute  for  W from  equation  6.6  to  get 


d i = - Fsr  sini  ( — gm- — ) cos  u d u . 

o g 


6.8 


6.8,  we  get  an  expression  for  d i. 


-119- 


d i = 


irR 


GM0(1 


u„R„  E R 

■]  t-2^)  1^)  (-2) 

e)  c mg  r 


x sine  sinicos  u d u 


6.9 


Since  i will  be  nearly  constant  over  one  orbit,  we  can  integrate  equa- 
tion 6.9  by  holding  sint  constant.  Also,  for  sine  we  must  substitute 
in  terms  of  u and  i.  The  change  in  inclination  after  one  orbit  is  then 


2 I ’ 2 2 

A i = - AiQ  ^ cos  u J 1 - sin  i sin  u d u , 


where 


o ttR  w R F R 

Ho  = T5  [ 2 ] ^ <s§>  <-?>  • 

0 ib  GMq(1  - 3)  C sd  r 


6.10 


s being  the  radius  of  the  particle  and  d its  density.  Also,  Fq  is  the 
energy  flux  of  solar  radiation  at  the  sun's  surface,  so  that 


Now  the  integral  in  u is  not  immediately  solvable,  but,  by 
expanding  the  term  under  the  square  root  into  a power  series  and 
integrating,  we  arrive  at 


Ai  = -ir  Af Q (1  - ^ sin^i  + -p;  sin^i  -...). 

Here  we  can  ignore  all  but  the  first  term  unless  i is  nearly  90°  (even 
then,  our  approximation  will  be  accurate  to  within  12%).  The  value 
for  Mq  can  be  found  by  substituting  the  appropriate  values  for  the 
constants  in  equation  6.10.  Notice  that  due  to  the  minus  sign  THIS  EFFECT 
ALWAYS  DECREASES  THE  INCLINATION.  For  example,  the  decrease  in  inclination 
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over  one  orbit,  for  s = 1 ym,  d = 1 -2^-,  and  sini  = is  0.28  seconds  of 
arc  per  orbit,  or  about  ^ minute  of  arc  for  the  100  orbits  a particle 
is  expected  to  survive  at  the  edge  of  the  evaporation  zone.  Obviously, 
this  mechanism  cannot  significantly  affect  the  inclination  of  dust 
particles  once  they  reach  the  F corona. 

We  can  compare  the  acceleration  on  the  particle  due  to  solar 
rotation  with  the  acceleration  on  a particle  near  the  sun  due  to  the 
strongest  perturber,  Jupiter  (see  Chapter  IV,  table  4.3).  The  ratio 
of  these  accelerations,  at  the  F corona,  is 

asr  _ sini 
aj  8sd  ’ 

with  d in  grams/cm  and  s in  microns.  A glance  at  table  4.1  shows 
that  this  effect  is  even  less  than  that  of  the  planets.  This  effect 
could  not  create  a density  enhancement  of  the  coronal  dust  in  the 
solar  equatorial  plane.  This  result  supports  our  conclusion  in  the 
analysis  of  the  observations,  performed  in  Chapter  V,  that  no  such 
enhancement  has  been  discovered  or  could  be  created  by  any  dynamical 
mechanism  in  the  F corona. 


CHAPTER  VII 


CONCLUSION 

The  F corona  is  a region  of  rapidly  varying  characteristics  in 
which  interplanetary  dust  is  subjected  to  numerous  dynamical  forces. 

We  have  attempted  to  clarify  the  effect  of  some  of  these  forces  on 

the  orbital  evolution  of  F coronal  dust,  using  both  theory  and  numerical 

simulations. 

In  Chapter  II  we  discovered  how  interplanetary  dust  distribution 
depends  on  solar  distance  through  interactions  with  the  inner  planets. 
Chapter  III  extended  the  distribution  of  interplanetary  dust  into  the 
F corona  as  a basis  for  the  initial  conditions  of  subsequent  simula- 
tions. In  Chapter  IV  we  derived  a method  for  determining  the  relative 
extent  of  each  planet's  zone  of  influence  when  perturbations  on  the 
orbits  of  particles  are  to  be  determined  and  used  this  to  show  that  the 
effect  of  planetary  perturbations  on  circumsolar  dust  is  negligible. 
Chapter  V was  an  in-depth  study  of  the  effect  of  the  solar  magnetic 
field  on  the  orbital  elements  of  dust  in  the  F corona.  Here,  we  have 
also  correlated  observations  of  the  F corona  by  previous  authors  with 
our  findings.  Finally  in  Chapter  VI  we  developed  an  expression  for 
the  radiation  pressure  force  of  a rotating  spherical  source. 

Among  our  findings  are  the  following: 

1.  Planetary  perturbations  increase  the  eccentricity  of  dust 

particle  orbits  and  decrease  the  P-R  lifetime  of  the  particles. 
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The  orbits  of  particles  which  do  not  experience  close  encounters 
with  the  planet  Jupiter  are,  on  the  average,  slowly  reduced  in 
inclination. 

2.  There  is  a strong  tendency  of  the  particles  to  align  their 
perihelia.  The  preferred  direction  is  toward  u>  = 20°. 

3.  Particles  enter  the  F corona  in  nearly  circular  orbits.  Although 
they  tend  toward  low  inclinations  due  to  planetary  perturba- 
tions, there  is  still  a large  spread  in  inclination  as  they 
reach  the  F corona.  We  found  no  preferred  direction  for  the 
longitudes  of  the  ascending  nodes  of  these  particles. 

4.  Planetary  perturbations  upon  F coronal  particles' orbits  are 
not  a significant  source  of  variation  in  the  orbital  elements. 

The  planet  Jupiter  is  the  strongest  perturber  in  the  F corona. 

5.  Interaction  with  the  solar  magnetic  field  does  not  affect  the 
semimajor  axis  or  eccentricity  of  dust  in  the  F corona.  It 
will  increase  the  inclinations  of  orbiting  dust  during  those 
times  when  the  equatorial  dipole  component  is  strong.  This 
effect  generates  a preferential  direction  along  the  dipole  axis 
in  the  line  of  nodes  of  the  orbiting  dust. 

6.  The  quadrupole  components  of  the  solar  magnetic  field  have  no 
effect  on  the  orbital  elements  of  circumsolar  dust. 

7.  Octupole  components  of  the  solar  magnetic  field  can  create 
equilibrium  positions  in  the  inclinations  and  the  longitudes  of 
the  ascending  nodes  of  dust  particle  orbits.  The  most  stable 
of  these  is  an  inclination  of  63°  from  the  g31  and  h31  terms. 

Only  the  g^  and  h^  octupole  terms  of  the  solar  magnetic  field 
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can  cause  an  alignment  of  the  interplanetary  dust  in  the  solar 
equatorial  plane,  and  this  effect  is  not  expected  to  play  a sig- 
nificant role  in  generating  features  in  the  solar  dust  cloud, 
since  the  magnitudes  of  the  octupole  components  rarely  exceed 
that  of  the  equatorial  dipole. 

8.  The  effect  of  solar  rotation  generates  a secular  decrease  in  the 
inclination  of  F coronal  dust,  but  this  change  is  very  small  and 
can  not  align  the  dust  into  a ring  structure. 

9.  Observations  of  the  F corona  imply  that  the  solar  dust  cloud 
extends  to  at  least  20°  above  and  below  the  solar  equatorial 
plane.  The  authors  who  contend  that  the  structure  is  a thin 
ring  have  yet  to  suggest  a mechanism  for  aligning  the  particle 
inclinations  into  a single  plane. 


CHAPTER  VIII 


SUGGESTIONS  FOR  FUTURE  RESEARCH 

Further  study  of  the  dynamical  mechanisms  of  the  F corona  may 
lead  to  a final  resolution  of  the  question  of  the  shape  of  the  solar 
dust  cloud.  Other  effects  to  study  include  the  following: 

1.  Alignment  of  particle  orbits  due  to  collisions  which  may  explain  the 
reduction  of  inclinations.  Collisions  also  play  an  important 

role  in  the  ultimate  lifetime  of  particles  in  the  F corona.  The 
lack  of  knowledge  of  the  exact  nature  of  these  particles  makes 
a meaningful  study  of  collision  processes  difficult  at  the 
present  time. 

2.  The  effect  of  sputtering  on  the  balance  between  P-R  drag  and  the 
sublimation  effect.  While  small  radial  and  tangential  forces 
such  as  those  due  to  the  magnetic  field  become  absorbed  in  this 
balance,  forces  due  to  sputtering  can  dominate  the  other  mech- 
anisms as  the  particle  passes  through  high  speed  streams. 

3.  Coulomb  drag  due  to  interaction  between  a charged  dust  grain  and 
the  solar  wind.  Like  the  previous  mechanism,  this  force  depends 
on  the  solar  wind  parameters  which  are  not  yet  well  understood 
in  this  region. 

4.  Anisotropies  in  radiation  pressure,  P-R  drag,  and  sputtering  due 
to  particle  rotation  or  surface  irregularities.  A study  of  the 
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effects  of  radiation  pressure  on  cylindrical  particles  (Voshchin- 
nikov  and  II 'in,  1983)  concluded  that  the  force  due  to  the  non- 
spherical  nature  of  the  particles  was  greater  than  P-R  drag.  A 
similar  result  is  obtained  from  the  Yarkovsky  effect  (the  aniso- 
tropic re-radiation  of  a rotating  particle;  Radzievskii,  1952). 
These  mechanisms  are  significant  even  outside  the  F corona. 

The  ultimate  answer  to  understanding  the  dynamics  of  dust  particles 
in  the  F corona  requires  long  period  observations  of  the  region  to 
determine  how  changes  in  the  structure  of  the  solar  dust  cloud  relate 
to  the  characteristics  of  the  solar  environment. 
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APPENDIX  A 


MATHEMATICAL  RELATIONS 


TRIGONOMETRIC  INTEGRALS 

A.  FIRST  ORDER  INTEGRALS 

<f  s i nu  d u 

B.  SECOND  ORDER  INTEGRALS 
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II.  LEGENDRE  POLYNOMIALS 

p®  = 1 
^0  1 

Pj  = cose 

P°  = \ (3cos2e-l) 

P3  = cose(5cos2e-3) 


pj  = sine 

1 2 2 
P2  = -3sinocose  P2  = 3sin 

13  2 2 
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P3  = - 1 5s i n3e 


P^  = (-l)mP7 
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l l 
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dmP^(cose) 

d(cose)m 


III.  COORDINATE  TRANSFORMATIONS 

Note:  The  angle  e is  the  polar  angle,  the  angle  $ is  the 

azimuthal  angle  as  measured  from  the  x-axis. 


x = sinecosifi  = cosn  cosu  - sinn  sinu  cosi 
y = sinesin<j>  = sinft  cosu  - cosfi  sinu  cosi 
z = cose  = sinu  sini 


IV.  GAUSSIAN  PERTURBATION  EQUATIONS 
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APPENDIX  B 


GRAPHIC  DATA  OF  PARTICLE  SIMULATIONS 
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Figure  B.l.  Number  density  vs.  inclination  at  0.35  A.U. 
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Figure  B.2.  Number  density  vs.  inclination  at  0.45  A.U. 
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Figure  B.4.  Number  density  vs.  inclination  at  0.65  A.U. 
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Figure  B.18.  Number  density  vs.  inclination  at 
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Figure  B.19.  Number  density  vs.  inclination  at  2.15  A.U. 


-153- 


<r 


H 


O 

►— < 


c 

LO 

CM 

CM 


C 

o 

-M 

fD 

c 


o 

c: 


LO 

> 


to 

c 

cd 

"O 

s- 

<D 

_Q 

E 

23 

Z 


CD 

CM 

CD 


CD 

S- 

3 

cn 


-154- 


<r 


o 


ZD 

c* 


LO 

O') 


OsJ 

4-> 

ro 


c 

o 


rO 


U 


CO 

> 

>> 

4-> 

•r- 

CO 

c 

CD 

"O 

i- 

<D 

JD 

E 

=3 


C\J 

CO 

CD 

13 

CD 


-155- 


-156- 


o 


LO 

LO 


CVJ 


4-> 

rO 


C 

o 


H3 


O 

c 


CO 

> 


to 

c 

a) 

XJ 

S~ 

aj 

jQ 


ro 

C\J 


CQ 


a; 

S- 


cn 


-157- 


LO 

UD 


C\J 

4-> 

ru 


C 

o 


03 

c 


o 

c 


CO 

> 

>> 

4- > 
•r— 

CO 

C 

CD 

-o 

5- 
Cl) 

-Q 

E 

13 


<3- 

CVJ 


CD 

S- 

13 

CT> 


-158- 


LO 

r- 


C\J 

+-> 

03 


c 

o 


03 

C 


o 

c 


to 

> 


to 

c 

CD 

X) 

S- 

CD 

_Q 

E 

3 


LO 

CNJ 

CD 

O) 

S- 

3 

cn 


100 


159- 


LO 

ro 


c 

o 


ro 

C 


u 

cz 


<U 

c 

LO 

o 

CJ 


(/) 

> 


LO 

C 

CD 

"O 


i- 

<D 

-O 


ViD 

C\J 


ca 


CD 

i- 

3 

CD 


100 


160 


Figure  B.27.  Number  density  vs.  cosine  of  inclination  at  0.45  A.U. 
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Figure  B.2S.  Number  density  vs.  cosine  of  inclination  at  0.65  A.U. 
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vs.  cosine  of  inclination  at  1.35  A.U. 
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Figure  B.39.  Number  density  vs.  cosine  of  inclination  at  1.65  A.U. 
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Figure  B.40.  Number  density  vs.  cosine  of  inclination  at  1.75  A.U. 
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Figure  B.41.  Number  density  vs.  cosine  of  inclination  at  1.85  A.U. 
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Figure  B.42.  Number  density  vs.  cosine  of  inclination  at  1.95  A.U. 
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Figure  B.45.  Number  density  vs.  cosine  of  inclination  at  2.25  A.U. 
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Figure  B.46.  Number  density  vs.  cosine  of  inclination  at  2.35  A.U. 
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Figure  B.50.  Number  density  vs.  cosine  of  inclination  at  2.75  A.U. 
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Figure  B.51.  Inclination  vs.  longi 
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Figure  B.52.  Inclination  vs.  longitude  of  the  ascending  node  at  0.45  A.U. 
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Figure  B.53.  Inclination  vs.  longitude  of  the  ascending  node  at  0.55  A.U. 
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Figure  B.54.  Inclination  vs.  longitude  of  the  ascending  node  at  0.65  A.U. 
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Figure  B.56.  Inclination  vs.  longitude  of  the  ascending  node  at  0.85  A.U. 
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Figure  B.57.  Inclination  vs.  longitude  of  the  ascending  node  at  0.95  A.U. 
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Figure  B.58.  Inclination  vs.  longitude  of  the  ascending  node  at  1.05  A.U. 
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igure  B.64.  Inclination  vs.  longitude  of  the  ascending  node  at  1.65  A.U. 
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Figure  B.73.  Inclination  vs.  longitude  of  the  ascending  node  at  2.55  A.U. 
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Figure  B.76.  Number  density  vs.  time  of  infall  at  0.35  A.U. 
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igure  B.81.  Number  density  vs.  time 
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Figure  B.85.  Number  density  vs.  time  of  infall  at  1.25  A.U. 
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Figure  B.86.  Number  density  vs.  time  of  infall  at  1.35  A.U. 
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Figure  B.89.  Number  density  vs.  time  of  infall  at  1.65  A.U. 
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igure  B.93.  Number  density  vs.  time 
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Figure  B.95.  Number  density  vs.  time  of  infall  at  2.25  A.U. 
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Figure  B.96.  Number  density  vs.  time  of  infall  at  2.35  A.U. 
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Figure  B.97.  Number  density  vs.  time  of  infali  at  2.45  A.U. 
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Figure  B.98.  Number  density  vs.  time 
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Figure  B.99.  Number  density  vs.  time  of  infail  at  2.65  A.U. 
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Figure  111.  Number  density  vs.  longitude  of  perihelion  at  1.35  A.U. 
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Figure  113.  Number  density  vs.  longitude  of  perihelion  at  1.55  A.U. 
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Figure  114.  Number  density  vs.  longitude  of  perihelion  at  1.65  A.U. 
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Figure  124.  Number  density  vs. 
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Figure  127.  Number  density  vs.  longitude  of  the  ascending  node  at  0.45  A.U. 
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Figure  131.  Number  density  vs.  longitude  of  the  ascending  node  at  0.85  A.U. 
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Figure  132.  Number  density  vs.  longitude  of  the  ascending  node  at  0. 
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Figure  133.  Number  density  vs.  longitude  of  the  ascending  node  at  1.05  A.U. 
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Figure  143.  Number  density  vs.  longitude  of  the  ascending  node  at  2.05  A.U. 
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Figure  149.  Number  density  vs.  longitude  of  the  ascending  node  at  2.65  A.U. 
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Figure  152.  Number  density  vs.  eccentricity  at  0.45  A.U. 
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Figure  155.  Number  density  vs.  eccentricity  at  0.75  A.U. 
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Figure  159.  Number  density  vs.  eccentricity  at  1.15  A.U. 


-293- 


>s> 


00  U?  ^ OJ 


o 


IXL 


o 

o 

UJ 


Figure  160.  Number  density  vs.  eccentricity  at  1.25  A.U. 
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igure  161.  Number  density  vs.  eccentricity  at  1.35  A.U. 
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Figure  167.  Number  density  vs.  eccentricity  at  1.95  A.U. 
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igure  168.  Number  density  vs.  eccentricity  at  2.05  A.U. 
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density  vs.  eccentricity  at  2.35  A.U. 
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Figure  174.  Number  density  vs.  eccentricity  at  2.65  A.U. 


10C.C 


308. 


-309- 


ce 

LU 

Q_ 


o 


o 

z 

o 


Figure  176.  Number  density  vs.  longitude  of  perihelion,  23,500  years  ago. 
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Figure  177.  Number  density  vs.  longitude  of  perihelion,  22,000  years  ago. 
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Figure  179.  Number  density  vs.  longitude  of  perihelion,  19,000  years  ago. 
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Figure  180.  Number  density  vs.  longitude  of  perihelion,  17,500  years  ago. 
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181.  Number  density  vs.  longitude  of  perihelion,  16,000  years  ago. 
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Figure  186.  Number  density  vs.  longitude  of  perihelion,  8,500  years  ago. 
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Figure  187.  Number  density  vs.  longitude  of  perihelion,  7,000  years  ago. 
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Figure  188.  Number  density  vs.  longitude  of  perihelion,  5,500  years  ago. 
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Figure  189.  Number  density  vs.  longitude  of  perihelion,  4,000  years  ago. 
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Figure  190.  Number  density  vs.  longitude  of  perihelion,  2,500  years  ago. 
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Figure  191.  Number  density  vs.  longitude  of  perihelion,  1,000  years  ago. 
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Figure  192.  Number  density  vs.  longitude  of  perihelion,  500  years  from  now. 
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Figure  193.  Number  density  vs.  longitude  of  perihelion,  2,000  years  from  now. 
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Figure  194.  Number  density  vs.  longitude  of  perihelion,  3,500  years  from  now. 
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